VIX Project : Forecasting Market Volatility¶

Table of Contents
¶

Overview
Importing Data, Functions and Packages
Exploring the Data

Addressing NaNs

How Measure Success
Modeling
1. ARCH/GARCH
2. VECM
3. AR
4. ARMA
Revised Strategy
Conclusion

Overview
¶

Data notes:¶

The data runs from March 2005 through end November 2016. It is daily.
The SPY is the S&P index and you have seen this one before.
The "RealVol" is realized volatility for the past 30 days i.e., the standard deviation of the returns on the S&P index in annualized form. See Dan's notebook on volatility.
The VIX is the CBOE volatility index. The current VIX index value quotes the expected annualized change in the S&P 500 index over the following 30 days, as computed from options-based theory and current options-market data. https://en.wikipedia.org/wiki/VIX
The VIXCM30 stands for the VIX 30 day volatility future, i.e., a measure of the VIX forward 30 days, i.e., from 30-60 days. Likewise the VIXCM60 is the 60-day forward future, i.e., a measure of the volatility between 60-90 days out in time. The five futures contracts constitute the "futures curve".

Project goals:¶

Imagine you are a trader and are interested in building a model for accurately forecasting the five futures prices, i.e., the five points in the volatility curve (VIXCM30, VIXCM60, VIXCM90, VIXCM120, VIXCM150). You can use all the data in the data frame above to come up with a forecasting model.
Devise a metric for this exercise. How would you convince someone with a good metric that the performance of your forecasting is highly accurate. There is no specific known metric, so feel free to look at existing research on volatility forecasting.
You may want to develop rolling forecasts, i.e., choose a reasonable look back period to fit the model and the a reasonable forward forecast. For example, you may want to use a 180-day look back to fit the model and then forecast the next weeks futures prices. (That is just an example, find your best approach.)

Importing Data, Functions and Packages
¶

First we will import the data and define functions which will be used throughout the notebook:¶

import pandas as pd
import matplotlib.pyplot as plt
from johansen_test import coint_johansen
from scipy import stats
%matplotlib inline
import numpy as np
from scipy import stats
from sklearn.model_selection import train_test_split
import statsmodels.api as sm # Time Series Analysis
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.metrics import mean_squared_error
import copy
import seaborn as sns
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 15,8
from statsmodels.tsa.ar_model import AR
from statsmodels.tsa.arima_model import ARIMA
import scipy
from scipy import *
import math
from scipy import interpolate
import arch
from arch import arch_model

from statsmodels.tsa.vector_ar.vecm import VECM
from sklearn.model_selection import train_test_split
import warnings
warnings.filterwarnings("ignore")

from statsmodels.tsa.stattools import adfuller

def run_ad_fuller(X):
    result = adfuller(X)
    print('ADF Statistic: %f' % result[0])
    print('p-value: %f' % result[1])
    print('Critical Values:')
    for key, value in result[4].items():
        print('\t%s: %.3f' % (key, value))
    if result[1] > 0.05:
        print('\nWe fail to reject the Null Hypothesis (H0) -- the time series has a unit root and is not stationary')
    else:
        print('\nWe reject the Null Hypothesis (H0) -- the time series is stationary')

from statsmodels.tsa.arima_model import ARMA
from prettytable import PrettyTable as pt
import time

# evaluate an ARMA model
def evaluate_arma_model(X, p, q, metric='BIC'):
    assert metric in ('AIC','BIC','HQIC','MSE')
    try:
        model = ARMA(X, order=(p,q))
        r = model.fit(trend='nc', method='css-mle')
        if metric == 'BIC':    value = r.bic
        elif metric == 'AIC':  value = r.aic
        elif metric == 'HQIC': value = r.hqic
        elif metric == 'MSE':  value = np.mean(r.sigma2) ** .5
        return value
    except:
        return None

def run_grid_over_arma(X, p_lags, q_lags, metric='BIC'):
    header = ['p\q'] + list(map(lambda q: 'q (%d)' % q, q_lags))
    start_time = time.time()
    table = pt(header)
    bad = []
    result = {}
    for p in p_lags:
        row = ['p(%d)' % p]
        for q in q_lags:
            #print(p,q)
            value = evaluate_arma_model(X, p, q, metric=metric)
            if value is None:
                bad += [(p,q)]
                row += ['n/a']
            else:
                row += ['%8.2F' % value]
                result[(p,q)] = value
                
        table.add_row(row)
        print(row)
    print(table)
    end_time = time.time()
    print('bad values: %s' % str(bad))
    print('time required', end_time - start_time)
    return result

from sklearn.metrics import mean_squared_error

def measure_error(actual, pred, label=None, show=True):
    sq_error = mean_squared_error(actual,pred)
    rmse = np.mean(sq_error)**.5
    avg = np.mean(actual)
    errors = {'RMSE': rmse, 'RMSE_pcent': rmse / avg, 'label':label}
    if show:
        for (k,v) in errors.items():
            print('%s = %s' % (k,v))
    return errors

def ts_plot(y, lags=None, figsize=(10, 8), style='bmh'):
    if not isinstance(y, pd.Series):
        y = pd.Series(y)
    with plt.style.context(style):    
        fig = plt.figure(figsize=figsize)
        layout = (3, 2)
        ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)
        acf_ax = plt.subplot2grid(layout, (1, 0))
        pacf_ax = plt.subplot2grid(layout, (1, 1))
        qq_ax = plt.subplot2grid(layout, (2, 0))
        pp_ax = plt.subplot2grid(layout, (2, 1))
        
        y.plot(ax=ts_ax)
        ts_ax.set_title('Time Series Analysis Plots')
        plot_acf(y, lags=lags, ax=acf_ax, alpha=0.5)
        plot_pacf(y, lags=lags, ax=pacf_ax, alpha=0.5)
        sm.qqplot(y, line='s', ax=qq_ax)
        qq_ax.set_title('QQ Plot')        
        stats.probplot(y, sparams=(y.mean(), y.std()), plot=pp_ax)

        plt.tight_layout()
    run_ad_fuller(y)
    return

from sklearn.metrics import confusion_matrix

def plot_confusion_matrix(actual, pred, show=True):
    y_true = []
    y_pred = []

    last = actual[0]
    for (a,p) in list(zip(actual, pred))[1:]:
        y_true += [-1 + int(a - last > 0) * 2]
        y_pred += [-1 + int(p - last > 0) * 2]
        last = a
    labels = ['tn', 'fp', 'fn', 'tp']
    cm = confusion_matrix(y_true, y_pred)
    
    print('Confusion Matrix\nRaw')
    print(cm)
    print('Normalized')
    cm_norm = cm / np.sum(cm)
    
    if show:
        print(cm_norm)
        plt.figure()
        plt.imshow(cm,  cmap=plt.cm.Blues)
        plt.colorbar()
        #tick_marks = np.arange(2)
        plt.xticks([-0.5,0.5], ['negative','positive'])
        plt.yticks([-0.5,0.5], ['negative','positive'])
        plt.ylabel('Actual')
        plt.xlabel('Predicted')
        plt.title('Manual AR(p) Model Predicted vs. Actual')
    return {'cm': cm, 'cm_norm': cm_norm}

from statsmodels.stats.diagnostic import acorr_ljungbox
def run_ljung_box(X,lags=40):
    (_,p_values) = acorr_ljungbox(X, lags=lags, boxpierce=False)
    print('Ljung-Box Test')
    for n,p_value in enumerate(p_values):
        if p_value < .05:
            print('lag: %3d | p_value: %.4f | still has information' % (n+1, p_value))

Importing the data¶

df_data=pd.read_csv('vix_data.csv')
df_data.head()

Exploring the Data
¶

Before we jump in we will make sure we understand the data available:¶

Volatility: The standard deviation of log returns where $ R = log(P_t / P_{t-1}) $
Realized Volatility: The volatility observer for an instrument over a recent window of N periods
- Note annualized daily volatility would be: $$ \sigma_{annual} = \sqrt{ 252 } * \sigma_{daily}$$
Futures: $F_0 = P_0$ + "cost of carry"
VIX: One of the most broadly used measure of volatility in the market
- The VIX calculation weighs the implied volatility of front month options (next monthly expiration to expire) and second month options (second month of monthly options to expire) option implied volatilities.
- The weights amount to a 30 day horizon of expiration.
VIX Futures: They represent forward expectations for volatility as well as the demand for insurance against tail events in the market
- In Contango: The VIX futures curve is upward sloping –> the spot VIX price is lower than the front month futures, which in turn are lower than the second month futures and so on.
- In Backwardation: The VIX futures curve is downward sloping– the spot VIX price is higher than the front month futures, which in turn are higher than the second month futures and so on.
- Note most of the time, VIX Futures are in Contango
Why are VIX Futures typically in Contango?:
- https://www.macroption.com/vix-futures-curve/
  - Contango can be interpreted in the way that the market expects the VIX index to increase from its current level going forward
  - Due to the skewed and mean reverting nature of the VIX and volatility in general (long time at low levels, with occasional big but mostly short-lived spikes) VIX Futures are typically in contango.
  - VIX futures contango is typically sharpest when the spot VIX Index is very low.
  - Backwardation typically occurs when the spot VIX index spikes and the market expects volatility to decrease again in the future.
- (http://www.bsam.com/volatility-futures-and-sp500-performance/)
  - The Black-Scholes model assumes that stock returns are distributed as a gaussian or normal distribution while actual market returns exhibit fat tails
  - Out of the money strike prices on both call and put options show higher implied volatility which reflects these fat tails
  - The actual frequency occurrence and magnitude of market crashes is larger than what a normal distribution would predict, which is now factored into option prices on the S&P500
  - The VIX futures theoretically provide a good asymmetric hedge against market risk since downside is limited and upside is substantial. This means that investors are likely willing to pay a premium to buy volatility futures to have some insurance against market drawdowns. That is why the volatility futures markets are typically in contango- reflecting this cost of insurance.

Below we can see how the VIX is mean reverting and skewed right (positively):

plt.rcParams['figure.figsize'] = 15,8
df_data.VIX.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x10a045be0>

df_data.VIX.hist(bins=100)

<matplotlib.axes._subplots.AxesSubplot at 0x1c19a4fdd8>

Now we will explore the data:¶

First let's summarize the data:

len(df_data)

2943

df_data.describe()

We can make a few observations from the summary statistics above:

VIXCM150 has many missing values, as there are only 509 out of 2943 rows filled
VIXCM120 also has some (but less) values missing
Looking at the mean values of VIX, VIXCM30, VIXCM60, VIXCM90, and VIXCM120 we can see that the VIX Futures seem to be in Contango. However, the VIXCM150 mean does not fall in line with this belief. This is possibly due to all of the missing values though.

Now we will also make some plots to visualize the data:

df_data.set_index('Date', inplace=True)
df_data.head()

Plotting all variables together:

df_data.VIXCM60.plot(label='VIXCM60')
df_data.VIX.plot(label='VIX')
df_data.VIXCM90.plot(label='VIXCM90')
df_data.VIXCM30.plot(label='VIXCM30')
df_data.VIXCM120.plot(label='VIXCM120')
df_data.VIXCM150.plot(label='VIXCM150')
df_data.RealVol.plot(label='RealVol')
df_data.SPY.plot(label='SPY')

<matplotlib.axes._subplots.AxesSubplot at 0x1c19c342b0>

Above we can see that SPY is upward trending and all of the VIX variables as well as the realized volatility may be cointegrated. We can take closer looks below:

df_data.SPY.plot(label='SPY')

<matplotlib.axes._subplots.AxesSubplot at 0x1c19c68c50>

df_data.VIXCM60.plot(label='VIXCM60')
df_data.VIX.plot(label='VIX')
df_data.VIXCM90.plot(label='VIXCM90')
df_data.VIXCM30.plot(label='VIXCM30')
df_data.VIXCM120.plot(label='VIXCM120')
df_data.VIXCM150.plot(label='VIXCM150')
df_data.RealVol.plot(label='RealVol')

<matplotlib.axes._subplots.AxesSubplot at 0x1c1c1e6668>

df_data.VIXCM60.plot(label='VIXCM60')
df_data.VIX.plot(label='VIX')
df_data.VIXCM90.plot(label='VIXCM90')
df_data.VIXCM30.plot(label='VIXCM30')
df_data.VIXCM120.plot(label='VIXCM120')
df_data.VIXCM150.plot(label='VIXCM150')

<matplotlib.axes._subplots.AxesSubplot at 0x1c1c423940>

We can see these value seem to be moving together and are likely cointegrated.

Addressing NaNs
¶

Addressing NAs¶

As we saw in our data exploration, both VIXCM120 and VIXCM150 have missing values.

VIXCM120, has only a handful of missing values, so we can look at where these values appear in the data:

df_data.isna().sum()

SPY            0
RealVol        0
VIX            0
VIXCM30        0
VIXCM60        0
VIXCM90        0
VIXCM120      28
VIXCM150    2434
dtype: int64

Since this variable is only missing data from before 3/8/2006, we could only use data after this date to avoid any issues the missing values could cause. But let's look at where the missing values in VIXCM150 are:

df_data['VIXCM150'].plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c1c6f27b8>

We can see these missing values are spread all over the data so we can't simply use a portion of the original dataset to avoid this issue- instead we will need to find another way to address the missing values for this variable. We were able to find a way to extrapolate points on a curve which we can use to fill in missing VIXCM120 and VIXCM150 values.

First we will use all of the VIX, VIXCM30, VIXCM60, and VIXCM90 values on a given day to extrapolate the VIXCM120 values on the missing days. Then we will use all of the VIX, VIXCM30, VIXCM60, VIXCM90 and VIXCM120 values on a given day to extrapolate the VIXCM150 values on the missing days.

We ran many tests to see how this extrapolation was performing compared to real VIXCM120 and VIXCM150 scores and we very happy this how accurate they were, so we believe this is a secure way to

for date,row in df_data.iterrows():
    if math.isnan(row['VIXCM120']):
        x=[0,30,60,90]
        y=[row['VIX'],row['VIXCM30'],row['VIXCM60'],row['VIXCM90']]
        f = interpolate.interp1d(x, y, fill_value = "extrapolate")
        new_y=float(f(120))
        df_data.set_value(date, 'VIXCM120', new_y)

for date,row in df_data.iterrows():
    if math.isnan(row['VIXCM150']):
        x=[0,30,60,90,120]
        y=[row['VIX'],row['VIXCM30'],row['VIXCM60'],row['VIXCM90'],row['VIXCM120']]
        f = interpolate.interp1d(x, y, fill_value = "extrapolate")
        new_y=float(f(150))
        df_data.set_value(date, 'VIXCM150', new_y)

df_data.isna().sum()

Date        0
SPY         0
RealVol     0
VIX         0
VIXCM30     0
VIXCM60     0
VIXCM90     0
VIXCM120    0
VIXCM150    0
dtype: int64

How Measure Success
¶

Before We Build Models We Will Establish How We Will Measure Their Success:¶

We will measure the success of a model with by observing the daily PnLs. For any given day, we will advise a direction, to buy (+1) or to sell (-1) a given future. Then we will multiply the direction by $Future Price_{t}$ - $Future Price_{t-1}$.

For example, if the VIXCM30 was 20 yesterday and the VIXCM30 is 24 today, and our model predicted the VIXCM30 would increase, then the PnL we'd output would be $(24-20)*(+1) = 4$ but if our model predicted the VIXCM30 would have decreased, our PnL output would have been $(24-20)*(-1) = -4$.

def make_money(actual,predicted):
    d = {'Predicted': predicted,
     'Actual': actual}
    PnL_df = pd.DataFrame(d)
    
    Total_PnL =0
    direction=0
    correct=0
    incorrect=0
    
    for i in range(0, (len(PnL_df)-1)):
        if PnL_df['Predicted'][i+1]>PnL_df['Actual'][i]:
            direction=1
        else:
            direction=-1
        daily_PnL= direction*(PnL_df.Actual[i+1]-PnL_df.Actual[i])
        Total_PnL+=daily_PnL
        if daily_PnL>=0:
            correct+=1
        else:
            incorrect+=1
            
    print('Total profit or loss: %.4f' % (Total_PnL))
    print('Number of correct days: %d' % (correct))
    print('total incorrect days: %d' % (incorrect))

Modeling
¶

Models To Be Tested:¶

1) ARCH/GARCH model

2) VECM Model

3) AR Models

4) ARMA Models

ARCH/GARCH
¶

ARCH and GARCH Models¶

With our ARCH and GARCH models we will try to predict the volatility of the SPY daily returns. Based on these predictions (volatility will increase or decrease), we can make a recommendation about whether to buy or sell each future.

First we will create a new variable for the SPY returns:

df_data['spy_ret'] = df_data['SPY'].pct_change() 
ts_plot(df_data['spy_ret'].dropna(), lags = 20)

ADF Statistic: -12.501615
p-value: 0.000000
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567

We reject the Null Hypothesis (H0) -- the time series is stationary

Next we can compute the variance of the returns. This variable will be used to see how accurate our model is:

# compute variance
spy_daily_avg = df_data['spy_ret'].mean()
df_data['variance'] = (df_data['spy_ret'] - spy_daily_avg) ** 2
ts_plot(df_data['variance'].dropna(), lags=20)

ADF Statistic: -4.766926
p-value: 0.000063
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567

We reject the Null Hypothesis (H0) -- the time series is stationary

Now we can run our model. Please note that we have run a number of ARCH and GARCH models (with varying orders) using the following code. We have summarized the results in a table after the code, which will just be run for one of the models.

df_data.reset_index(inplace=True)

train,test = train_test_split(df_data, test_size=0.2, shuffle=False)

df_data['Date'] = pd.to_datetime(df_data['Date'])
lst = [0,29,59,89,119,149] # VIXCM nth values

H = train['spy_ret'].dropna()[-180:].tolist()
#H = train['spy_logret'].dropna()[-180:].tolist()
ACTUAL = []
PRED = []
for date,row in test.iterrows():
    am = arch_model(H, vol='Arch', p=1, mean='Constant')
    #am = arch_model(H, vol='Garch', p=3,q=2, dist='Normal',mean='Constant')
    f = am.fit()
    fc = f.forecast(horizon=150)
    ACTUAL += [row['variance']]
    PRED += [fc.variance.values[-1][lst]]
    #print(date,ACTUAL[-1],PRED[-1])
    H += [row['spy_ret']]
    H = H[1:]

Iteration:      1,   Func. Count:      5,   Neg. LLF: -659.9420976322722
Iteration:      2,   Func. Count:     17,   Neg. LLF: -660.0753368434889
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -660.0753377375784
            Iterations: 3
            Function evaluations: 25
            Gradient evaluations: 2
Iteration:      1,   Func. Count:      5,   Neg. LLF: -662.3038993947875
Iteration:      2,   Func. Count:     17,   Neg. LLF: -662.4368207461009
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -662.4368214190276
            Iterations: 3
            Function evaluations: 25
            Gradient evaluations: 2
Iteration:      1,   Func. Count:      5,   Neg. LLF: -655.646871393028
Iteration:      2,   Func. Count:     17,   Neg. LLF: -655.7033109942909
Iteration:      3,   Func. Count:     24,   Neg. LLF: -655.7286575610116
Iteration:      4,   Func. Count:     33,   Neg. LLF: -655.7383757791745
Iteration:      5,   Func. Count:     45,   Neg. LLF: -655.7383822997087
Iteration:      6,   Func. Count:     53,   Neg. LLF: -655.7397939858729
Iteration:      7,   Func. Count:     63,   Neg. LLF: -655.7398118266389
Iteration:      8,   Func. Count:     71,   Neg. LLF: -655.7555757455966
Iteration:      9,   Func. Count:     78,   Neg. LLF: -655.7618562021332
Iteration:     10,   Func. Count:     85,   Neg. LLF: -655.7780134636225
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -655.7780135141056
            Iterations: 11
            Function evaluations: 91
            Gradient evaluations: 10
Iteration:      1,   Func. Count:      5,   Neg. LLF: -655.2328342761537
Iteration:      2,   Func. Count:     17,   Neg. LLF: -655.2901630573631
Iteration:      3,   Func. Count:     25,   Neg. LLF: -655.2913706138233
Iteration:      4,   Func. Count:     33,   Neg. LLF: -655.2943181243002
Iteration:      5,   Func. Count:     42,   Neg. LLF: -655.2977253321736
Iteration:      6,   Func. Count:     49,   Neg. LLF: -655.311452941711
Iteration:      7,   Func. Count:     56,   Neg. LLF: -655.3114718372431
Iteration:      8,   Func. Count:     61,   Neg. LLF: -655.3115187490898
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -655.3115187491103
            Iterations: 8
            Function evaluations: 61
            Gradient evaluations: 8
Iteration:      1,   Func. Count:      5,   Neg. LLF: -655.3916172522531
Iteration:      2,   Func. Count:     17,   Neg. LLF: -655.450785586859
Iteration:      3,   Func. Count:     25,   Neg. LLF: -655.4510087778469
Iteration:      4,   Func. Count:     33,   Neg. LLF: -655.4562012432839
Iteration:      5,   Func. Count:     40,   Neg. LLF: -655.4750025166768
Iteration:      6,   Func. Count:     48,   Neg. LLF: -655.475054695914
Iteration:      7,   Func. Count:     53,   Neg. LLF: -655.4751238131926
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -655.4751243290468
            Iterations: 7
            Function evaluations: 55
            Gradient evaluations: 7
Iteration:      1,   Func. Count:      5,   Neg. LLF: -654.3000340890056
Iteration:      2,   Func. Count:     17,   Neg. LLF: -654.3551449772781
Iteration:      3,   Func. Count:     25,   Neg. LLF: -654.3622395843736
Iteration:      4,   Func. Count:     32,   Neg. LLF: -654.3802359039512
Iteration:      5,   Func. Count:     39,   Neg. LLF: -654.3802382805384
Iteration:      6,   Func. Count:     46,   Neg. LLF: -654.380270876936
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -654.3802708769763
            Iterations: 6
            Function evaluations: 46
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -654.330742040216
Iteration:      2,   Func. Count:     17,   Neg. LLF: -654.382317278534
Iteration:      3,   Func. Count:     26,   Neg. LLF: -654.3823288895812
Iteration:      4,   Func. Count:     34,   Neg. LLF: -654.3873314189232
Iteration:      5,   Func. Count:     43,   Neg. LLF: -654.3895391000727
Iteration:      6,   Func. Count:     50,   Neg. LLF: -654.3914319022268
Iteration:      7,   Func. Count:     57,   Neg. LLF: -654.3914708832925
Iteration:      8,   Func. Count:     62,   Neg. LLF: -654.3914956555268
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -654.3914956555341
            Iterations: 8
            Function evaluations: 62
            Gradient evaluations: 8
Iteration:      1,   Func. Count:      5,   Neg. LLF: -654.0680900349303
Iteration:      2,   Func. Count:     17,   Neg. LLF: -654.1194724827947
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -654.1194724653178
            Iterations: 3
            Function evaluations: 26
            Gradient evaluations: 2
Iteration:      1,   Func. Count:      5,   Neg. LLF: -652.6787825517026
Iteration:      2,   Func. Count:     17,   Neg. LLF: -652.7336654588219
Iteration:      3,   Func. Count:     27,   Neg. LLF: -652.7336673825902
Iteration:      4,   Func. Count:     39,   Neg. LLF: -652.7336731675352
Iteration:      5,   Func. Count:     47,   Neg. LLF: -652.7338865653126
Iteration:      6,   Func. Count:     57,   Neg. LLF: -652.7338936860979
Iteration:      7,   Func. Count:     69,   Neg. LLF: -652.733944285272
Iteration:      8,   Func. Count:     77,   Neg. LLF: -652.7340009741038
Iteration:      9,   Func. Count:     87,   Neg. LLF: -652.734181995491
Iteration:     10,   Func. Count:     95,   Neg. LLF: -652.7446197395984
Iteration:     11,   Func. Count:    102,   Neg. LLF: -652.7470880450417
Iteration:     12,   Func. Count:    110,   Neg. LLF: -652.7472135420189
Iteration:     13,   Func. Count:    116,   Neg. LLF: -652.747314124274
Iteration:     14,   Func. Count:    121,   Neg. LLF: -652.7473408847382
Iteration:     15,   Func. Count:    126,   Neg. LLF: -652.7473482494086
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -652.7473482494088
            Iterations: 17
            Function evaluations: 126
            Gradient evaluations: 15
Iteration:      1,   Func. Count:      5,   Neg. LLF: -652.6926809215587
Iteration:      2,   Func. Count:     17,   Neg. LLF: -652.7449078342412
Iteration:      3,   Func. Count:     25,   Neg. LLF: -652.751888295217
Iteration:      4,   Func. Count:     34,   Neg. LLF: -652.7535859687284
Iteration:      5,   Func. Count:     42,   Neg. LLF: -652.753762089377
Iteration:      6,   Func. Count:     50,   Neg. LLF: -652.7538116508761
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -652.753811650862
            Iterations: 7
            Function evaluations: 61
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -653.3717631080954
Iteration:      2,   Func. Count:     17,   Neg. LLF: -653.4268053312135
Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: -653.4268053133221
            Iterations: 6
            Function evaluations: 17
            Gradient evaluations: 2
Iteration:      1,   Func. Count:      5,   Neg. LLF: -653.1205996925453
Iteration:      2,   Func. Count:     17,   Neg. LLF: -653.1756980444384
Iteration:      3,   Func. Count:     25,   Neg. LLF: -653.185206935576
Iteration:      4,   Func. Count:     33,   Neg. LLF: -653.1855298948992
Iteration:      5,   Func. Count:     41,   Neg. LLF: -653.1855397252682
Iteration:      6,   Func. Count:     46,   Neg. LLF: -653.1855443275228
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -653.1855443275244
            Iterations: 6
            Function evaluations: 46
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -652.8880277736112
Iteration:      2,   Func. Count:     17,   Neg. LLF: -652.9443040753004
Iteration:      3,   Func. Count:     25,   Neg. LLF: -652.9512913092688
Iteration:      4,   Func. Count:     34,   Neg. LLF: -652.9539662079872
Iteration:      5,   Func. Count:     41,   Neg. LLF: -652.95400031495
Iteration:      6,   Func. Count:     49,   Neg. LLF: -652.9540516246801
Iteration:      7,   Func. Count:     55,   Neg. LLF: -652.9540545809866
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -652.9540545811503
            Iterations: 7
            Function evaluations: 66
            Gradient evaluations: 7

/Users/austincostello1/anaconda3/lib/python3.6/site-packages/arch/univariate/base.py:522: ConvergenceWarning: 
The optimizer returned code 8. The message is:
Positive directional derivative for linesearch
See scipy.optimize.fmin_slsqp for code meaning.

  ConvergenceWarning)

Iteration:      1,   Func. Count:      5,   Neg. LLF: -652.8470560679114
Iteration:      2,   Func. Count:     17,   Neg. LLF: -652.9047980862209
Iteration:      3,   Func. Count:     25,   Neg. LLF: -652.9132030394405
Iteration:      4,   Func. Count:     34,   Neg. LLF: -652.9147329200343
Iteration:      5,   Func. Count:     41,   Neg. LLF: -652.915463318544
Iteration:      6,   Func. Count:     49,   Neg. LLF: -652.9155194860964
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -652.9155194860865
            Iterations: 6
            Function evaluations: 49
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -652.0762467825074
Iteration:      2,   Func. Count:     17,   Neg. LLF: -652.1071195416826
Iteration:      3,   Func. Count:     26,   Neg. LLF: -652.1076621772175
Iteration:      4,   Func. Count:     33,   Neg. LLF: -652.1526014125025
Iteration:      5,   Func. Count:     41,   Neg. LLF: -652.1527779284752
Iteration:      6,   Func. Count:     46,   Neg. LLF: -652.1529879321919
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -652.1529879322479
            Iterations: 6
            Function evaluations: 46
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -651.8083581854055
Iteration:      2,   Func. Count:     17,   Neg. LLF: -651.8400072764684
Iteration:      3,   Func. Count:     26,   Neg. LLF: -651.8404942369202
Iteration:      4,   Func. Count:     33,   Neg. LLF: -651.874945906507
Iteration:      5,   Func. Count:     42,   Neg. LLF: -651.8751517410087
Iteration:      6,   Func. Count:     47,   Neg. LLF: -651.8752237767426
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -651.8752237768247
            Iterations: 6
            Function evaluations: 47
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -651.9045345174331
Iteration:      2,   Func. Count:     17,   Neg. LLF: -651.9385393875501
Iteration:      3,   Func. Count:     26,   Neg. LLF: -651.9388898606351
Iteration:      4,   Func. Count:     33,   Neg. LLF: -651.9806663293878
Iteration:      5,   Func. Count:     41,   Neg. LLF: -651.9808033301526
Iteration:      6,   Func. Count:     46,   Neg. LLF: -651.981041529755
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -651.9810415298115
            Iterations: 6
            Function evaluations: 46
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -652.1866631561422
Iteration:      2,   Func. Count:     17,   Neg. LLF: -652.2553018455044
Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: -652.255301838319
            Iterations: 6
            Function evaluations: 17
            Gradient evaluations: 2
Iteration:      1,   Func. Count:      5,   Neg. LLF: -652.1668330825432
Iteration:      2,   Func. Count:     17,   Neg. LLF: -652.2401795888227
Iteration:      3,   Func. Count:     25,   Neg. LLF: -652.2470896837203
Iteration:      4,   Func. Count:     34,   Neg. LLF: -652.2499254369105
Iteration:      5,   Func. Count:     41,   Neg. LLF: -652.2503234580743
Iteration:      6,   Func. Count:     49,   Neg. LLF: -652.2503705854149
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -652.2503705854194
            Iterations: 6
            Function evaluations: 49
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -652.4054519966676
Iteration:      2,   Func. Count:     17,   Neg. LLF: -652.4668739241156
Iteration:      3,   Func. Count:     26,   Neg. LLF: -652.4668993117391
Iteration:      4,   Func. Count:     35,   Neg. LLF: -652.4687251838101
Iteration:      5,   Func. Count:     43,   Neg. LLF: -652.4726293166317
Iteration:      6,   Func. Count:     52,   Neg. LLF: -652.4732334398202
Iteration:      7,   Func. Count:     59,   Neg. LLF: -652.4757031138633
Iteration:      8,   Func. Count:     67,   Neg. LLF: -652.476840989589
Iteration:      9,   Func. Count:     72,   Neg. LLF: -652.4768781785262
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -652.4768781785278
            Iterations: 9
            Function evaluations: 72
            Gradient evaluations: 9

/Users/austincostello1/anaconda3/lib/python3.6/site-packages/arch/univariate/base.py:522: ConvergenceWarning: 
The optimizer returned code 8. The message is:
Positive directional derivative for linesearch
See scipy.optimize.fmin_slsqp for code meaning.

  ConvergenceWarning)

Iteration:      1,   Func. Count:      5,   Neg. LLF: -653.7659050297555
Iteration:      2,   Func. Count:     17,   Neg. LLF: -653.8339200601217
Iteration:      3,   Func. Count:     25,   Neg. LLF: -653.8340661641085
Iteration:      4,   Func. Count:     33,   Neg. LLF: -653.8368287400006
Iteration:      5,   Func. Count:     41,   Neg. LLF: -653.8417374193965
Iteration:      6,   Func. Count:     48,   Neg. LLF: -653.8423238619855
Iteration:      7,   Func. Count:     56,   Neg. LLF: -653.8430647547218
Iteration:      8,   Func. Count:     61,   Neg. LLF: -653.8434126869214
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -653.8434126868774
            Iterations: 8
            Function evaluations: 61
            Gradient evaluations: 8
Iteration:      1,   Func. Count:      5,   Neg. LLF: -653.8218614908806
Iteration:      2,   Func. Count:     17,   Neg. LLF: -653.893208460891
Iteration:      3,   Func. Count:     26,   Neg. LLF: -653.8933163918327
Iteration:      4,   Func. Count:     34,   Neg. LLF: -653.8968445258797
Iteration:      5,   Func. Count:     43,   Neg. LLF: -653.8993177530855
Iteration:      6,   Func. Count:     50,   Neg. LLF: -653.9049323382559
Iteration:      7,   Func. Count:     58,   Neg. LLF: -653.9049356127038
Iteration:      8,   Func. Count:     63,   Neg. LLF: -653.904950670782
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -653.9049506707827
            Iterations: 8
            Function evaluations: 63
            Gradient evaluations: 8
Iteration:      1,   Func. Count:      5,   Neg. LLF: -654.2775031840852
Iteration:      2,   Func. Count:     17,   Neg. LLF: -654.3438191817922
Iteration:      3,   Func. Count:     24,   Neg. LLF: -654.3748976068796
Iteration:      4,   Func. Count:     32,   Neg. LLF: -654.382716608664
Iteration:      5,   Func. Count:     40,   Neg. LLF: -654.3894044831445
Iteration:      6,   Func. Count:     47,   Neg. LLF: -654.3923874351515
Iteration:      7,   Func. Count:     55,   Neg. LLF: -654.3925100754279
Iteration:      8,   Func. Count:     60,   Neg. LLF: -654.3925340873183
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -654.3925345583884
            Iterations: 8
            Function evaluations: 62
            Gradient evaluations: 8
Iteration:      1,   Func. Count:      5,   Neg. LLF: -656.1309677044703
Iteration:      2,   Func. Count:     17,   Neg. LLF: -656.2240946360441
Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: -656.2240946184189
            Iterations: 6
            Function evaluations: 17
            Gradient evaluations: 2
Iteration:      1,   Func. Count:      5,   Neg. LLF: -656.3475186553014
Iteration:      2,   Func. Count:     17,   Neg. LLF: -656.4439621790955
Iteration:      3,   Func. Count:     26,   Neg. LLF: -656.4443720977329
Iteration:      4,   Func. Count:     35,   Neg. LLF: -656.4486198996442
Iteration:      5,   Func. Count:     44,   Neg. LLF: -656.4493039245647
Iteration:      6,   Func. Count:     56,   Neg. LLF: -656.4510582531769
Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: -656.4510582368882
            Iterations: 10
            Function evaluations: 56
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -656.4275507400528
Iteration:      2,   Func. Count:     17,   Neg. LLF: -656.5140017245888
Iteration:      3,   Func. Count:     26,   Neg. LLF: -656.5185967945338
Iteration:      4,   Func. Count:     34,   Neg. LLF: -656.5317035188066
Iteration:      5,   Func. Count:     42,   Neg. LLF: -656.537977696685
Iteration:      6,   Func. Count:     51,   Neg. LLF: -656.5380521784421
Iteration:      7,   Func. Count:     58,   Neg. LLF: -656.5384777179393
Iteration:      8,   Func. Count:     66,   Neg. LLF: -656.5388948878749
Iteration:      9,   Func. Count:     71,   Neg. LLF: -656.5391280236668
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -656.5391280236361
            Iterations: 10
            Function evaluations: 71
            Gradient evaluations: 9
Iteration:      1,   Func. Count:      5,   Neg. LLF: -656.7588416212634
Iteration:      2,   Func. Count:     17,   Neg. LLF: -656.8483431955644
Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: -656.8483431711688
            Iterations: 6
            Function evaluations: 17
            Gradient evaluations: 2

/Users/austincostello1/anaconda3/lib/python3.6/site-packages/arch/univariate/base.py:522: ConvergenceWarning: 
The optimizer returned code 8. The message is:
Positive directional derivative for linesearch
See scipy.optimize.fmin_slsqp for code meaning.

  ConvergenceWarning)
/Users/austincostello1/anaconda3/lib/python3.6/site-packages/arch/univariate/base.py:522: ConvergenceWarning: 
The optimizer returned code 8. The message is:
Positive directional derivative for linesearch
See scipy.optimize.fmin_slsqp for code meaning.

  ConvergenceWarning)
/Users/austincostello1/anaconda3/lib/python3.6/site-packages/arch/univariate/base.py:522: ConvergenceWarning: 
The optimizer returned code 8. The message is:
Positive directional derivative for linesearch
See scipy.optimize.fmin_slsqp for code meaning.

  ConvergenceWarning)

Iteration:      1,   Func. Count:      5,   Neg. LLF: -657.5005382225817
Iteration:      2,   Func. Count:     17,   Neg. LLF: -657.5631829228405
Iteration:      3,   Func. Count:     25,   Neg. LLF: -657.6160720755558
Iteration:      4,   Func. Count:     33,   Neg. LLF: -657.6194177486296
Iteration:      5,   Func. Count:     40,   Neg. LLF: -657.6202661269544
Iteration:      6,   Func. Count:     48,   Neg. LLF: -657.6212855551237
Iteration:      7,   Func. Count:     53,   Neg. LLF: -657.6213314992976
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -657.6213320053339
            Iterations: 7
            Function evaluations: 55
            Gradient evaluations: 7
Iteration:      1,   Func. Count:      5,   Neg. LLF: -660.8030295819241
Iteration:      2,   Func. Count:     17,   Neg. LLF: -660.8612062282516
Iteration:      3,   Func. Count:     25,   Neg. LLF: -660.861611824671
Iteration:      4,   Func. Count:     33,   Neg. LLF: -660.8667693345633
Iteration:      5,   Func. Count:     42,   Neg. LLF: -660.8693858824538
Iteration:      6,   Func. Count:     49,   Neg. LLF: -660.9072184098063
Iteration:      7,   Func. Count:     58,   Neg. LLF: -660.9072449966604
Iteration:      8,   Func. Count:     63,   Neg. LLF: -660.9072733124492
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -660.9072733124756
            Iterations: 8
            Function evaluations: 63
            Gradient evaluations: 8
Iteration:      1,   Func. Count:      5,   Neg. LLF: -660.1572169625075
Iteration:      2,   Func. Count:     17,   Neg. LLF: -660.2106211711047
Iteration:      3,   Func. Count:     25,   Neg. LLF: -660.210988626591
Iteration:      4,   Func. Count:     33,   Neg. LLF: -660.216770861872
Iteration:      5,   Func. Count:     42,   Neg. LLF: -660.2203082867882
Iteration:      6,   Func. Count:     49,   Neg. LLF: -660.255652517247
Iteration:      7,   Func. Count:     57,   Neg. LLF: -660.2556607954765
Iteration:      8,   Func. Count:     62,   Neg. LLF: -660.2556975852559
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -660.2556975852802
            Iterations: 8
            Function evaluations: 62
            Gradient evaluations: 8
Iteration:      1,   Func. Count:      5,   Neg. LLF: -660.3320532461248
Iteration:      2,   Func. Count:     17,   Neg. LLF: -660.3798514557893
Iteration:      3,   Func. Count:     25,   Neg. LLF: -660.3883133013835
Iteration:      4,   Func. Count:     32,   Neg. LLF: -660.4141464890415
Iteration:      5,   Func. Count:     41,   Neg. LLF: -660.4141476923797
Iteration:      6,   Func. Count:     46,   Neg. LLF: -660.4141825416148
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -660.4141825416422
            Iterations: 6
            Function evaluations: 46
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -660.6231341153925
Iteration:      2,   Func. Count:     17,   Neg. LLF: -660.6727945035133
Iteration:      3,   Func. Count:     25,   Neg. LLF: -660.680662287278
Iteration:      4,   Func. Count:     32,   Neg. LLF: -660.6894251618182
Iteration:      5,   Func. Count:     39,   Neg. LLF: -660.7063410006107
Iteration:      6,   Func. Count:     47,   Neg. LLF: -660.7064102343022
Iteration:      7,   Func. Count:     52,   Neg. LLF: -660.7064132307494
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -660.7064132306159
            Iterations: 7
            Function evaluations: 63
            Gradient evaluations: 7
Iteration:      1,   Func. Count:      5,   Neg. LLF: -660.0493624855969
Iteration:      2,   Func. Count:     17,   Neg. LLF: -660.0931453024308
Iteration:      3,   Func. Count:     26,   Neg. LLF: -660.0932014960817
Iteration:      4,   Func. Count:     38,   Neg. LLF: -660.0932601130977
Iteration:      5,   Func. Count:     45,   Neg. LLF: -660.0983403396945
Iteration:      6,   Func. Count:     55,   Neg. LLF: -660.0985312738595
Iteration:      7,   Func. Count:     64,   Neg. LLF: -660.1032107134483
Iteration:      8,   Func. Count:     76,   Neg. LLF: -660.106772968092
Iteration:      9,   Func. Count:     84,   Neg. LLF: -660.1163705402789
Iteration:     10,   Func. Count:     91,   Neg. LLF: -660.1205298917457
Iteration:     11,   Func. Count:     98,   Neg. LLF: -660.1220507108571
Iteration:     12,   Func. Count:    103,   Neg. LLF: -660.1229598646087
Iteration:     13,   Func. Count:    110,   Neg. LLF: -660.1229622918388
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -660.1229622917242
            Iterations: 15
            Function evaluations: 110
            Gradient evaluations: 13
Iteration:      1,   Func. Count:      5,   Neg. LLF: -660.1985093070596
Iteration:      2,   Func. Count:     17,   Neg. LLF: -660.2395942810831
Iteration:      3,   Func. Count:     24,   Neg. LLF: -660.2400789573148
Iteration:      4,   Func. Count:     36,   Neg. LLF: -660.2401139372516
Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: -660.2401139134072
            Iterations: 8
            Function evaluations: 36
            Gradient evaluations: 4
Iteration:      1,   Func. Count:      5,   Neg. LLF: -659.4560195330346
Iteration:      2,   Func. Count:     17,   Neg. LLF: -659.4997416573165
Iteration:      3,   Func. Count:     26,   Neg. LLF: -659.5009224971657
Iteration:      4,   Func. Count:     35,   Neg. LLF: -659.504412086683
Iteration:      5,   Func. Count:     44,   Neg. LLF: -659.5079353134158
Iteration:      6,   Func. Count:     56,   Neg. LLF: -659.5088849114567
Iteration:      7,   Func. Count:     64,   Neg. LLF: -659.5089982893812
Iteration:      8,   Func. Count:     73,   Neg. LLF: -659.5106710282314
Iteration:      9,   Func. Count:     80,   Neg. LLF: -659.5139590356489
Iteration:     10,   Func. Count:     87,   Neg. LLF: -659.5146103200287
Iteration:     11,   Func. Count:     94,   Neg. LLF: -659.5182142979273
Iteration:     12,   Func. Count:    101,   Neg. LLF: -659.5184419400191
Iteration:     13,   Func. Count:    106,   Neg. LLF: -659.5186368903434
Iteration:     14,   Func. Count:    111,   Neg. LLF: -659.5186735433208
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -659.518673543316
            Iterations: 16
            Function evaluations: 111
            Gradient evaluations: 14
Iteration:      1,   Func. Count:      5,   Neg. LLF: -659.3930773033558
Iteration:      2,   Func. Count:     17,   Neg. LLF: -659.435800469605
Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: -659.435800449866
            Iterations: 6
            Function evaluations: 17
            Gradient evaluations: 2

/Users/austincostello1/anaconda3/lib/python3.6/site-packages/arch/univariate/base.py:522: ConvergenceWarning: 
The optimizer returned code 8. The message is:
Positive directional derivative for linesearch
See scipy.optimize.fmin_slsqp for code meaning.

  ConvergenceWarning)
/Users/austincostello1/anaconda3/lib/python3.6/site-packages/arch/univariate/base.py:522: ConvergenceWarning: 
The optimizer returned code 8. The message is:
Positive directional derivative for linesearch
See scipy.optimize.fmin_slsqp for code meaning.

  ConvergenceWarning)

Iteration:      1,   Func. Count:      5,   Neg. LLF: -659.5776210689435
Iteration:      2,   Func. Count:     17,   Neg. LLF: -659.6405081058016
Iteration:      3,   Func. Count:     25,   Neg. LLF: -659.6498544451362
Iteration:      4,   Func. Count:     32,   Neg. LLF: -659.6668815752219
Iteration:      5,   Func. Count:     40,   Neg. LLF: -659.6673872894473
Iteration:      6,   Func. Count:     46,   Neg. LLF: -659.6679996187197
Iteration:      7,   Func. Count:     51,   Neg. LLF: -659.6680040415179
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -659.6680040415159
            Iterations: 7
            Function evaluations: 51
            Gradient evaluations: 7
Iteration:      1,   Func. Count:      5,   Neg. LLF: -660.94979376385
Iteration:      2,   Func. Count:     17,   Neg. LLF: -661.0129741055507
Iteration:      3,   Func. Count:     25,   Neg. LLF: -661.0216044556266
Iteration:      4,   Func. Count:     32,   Neg. LLF: -661.0262543460963
Iteration:      5,   Func. Count:     39,   Neg. LLF: -661.0361382103742
Iteration:      6,   Func. Count:     48,   Neg. LLF: -661.0361421413608
Iteration:      7,   Func. Count:     53,   Neg. LLF: -661.0361575095683
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -661.0361575095288
            Iterations: 7
            Function evaluations: 53
            Gradient evaluations: 7
Iteration:      1,   Func. Count:      5,   Neg. LLF: -659.8837409512778
Iteration:      2,   Func. Count:     17,   Neg. LLF: -659.9502072252587
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -659.9502072897734
            Iterations: 2
            Function evaluations: 24
            Gradient evaluations: 2
Iteration:      1,   Func. Count:      5,   Neg. LLF: -659.5339019562297
Iteration:      2,   Func. Count:     17,   Neg. LLF: -659.5973975488157
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -659.5973984699548
            Iterations: 2
            Function evaluations: 23
            Gradient evaluations: 2
Iteration:      1,   Func. Count:      5,   Neg. LLF: -659.1988066725972
Iteration:      2,   Func. Count:     17,   Neg. LLF: -659.2607100999071
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -659.2607100753769
            Iterations: 2
            Function evaluations: 26
            Gradient evaluations: 2
Iteration:      1,   Func. Count:      5,   Neg. LLF: -655.8021559108395
Iteration:      2,   Func. Count:     17,   Neg. LLF: -655.8472654895512
Iteration:      3,   Func. Count:     25,   Neg. LLF: -655.8507894489
Iteration:      4,   Func. Count:     37,   Neg. LLF: -655.8603993662937
Iteration:      5,   Func. Count:     45,   Neg. LLF: -655.8661742425697
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -655.8661743672712
            Iterations: 7
            Function evaluations: 53
            Gradient evaluations: 5
Iteration:      1,   Func. Count:      5,   Neg. LLF: -655.0436978709222
Iteration:      2,   Func. Count:     17,   Neg. LLF: -655.0696426843193
Iteration:      3,   Func. Count:     26,   Neg. LLF: -655.0697052161916
Iteration:      4,   Func. Count:     35,   Neg. LLF: -655.0711767067841
Iteration:      5,   Func. Count:     43,   Neg. LLF: -655.0828128433426
Iteration:      6,   Func. Count:     50,   Neg. LLF: -655.0911165812582
Iteration:      7,   Func. Count:     58,   Neg. LLF: -655.0911328286559
Iteration:      8,   Func. Count:     64,   Neg. LLF: -655.0925077051585
Iteration:      9,   Func. Count:     69,   Neg. LLF: -655.0925104495407
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -655.0925104495443
            Iterations: 9
            Function evaluations: 69
            Gradient evaluations: 9
Iteration:      1,   Func. Count:      5,   Neg. LLF: -655.8378448819667
Iteration:      2,   Func. Count:     17,   Neg. LLF: -655.9219989316725
Iteration:      3,   Func. Count:     26,   Neg. LLF: -655.9220414887947
Iteration:      4,   Func. Count:     35,   Neg. LLF: -655.9228590820528
Iteration:      5,   Func. Count:     44,   Neg. LLF: -655.925783667435
Iteration:      6,   Func. Count:     56,   Neg. LLF: -655.9264987625086
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -655.9264996647587
            Iterations: 7
            Function evaluations: 62
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -657.6172047154266
Iteration:      2,   Func. Count:     17,   Neg. LLF: -657.6876351008295
Iteration:      3,   Func. Count:     26,   Neg. LLF: -657.6902978981756
Iteration:      4,   Func. Count:     38,   Neg. LLF: -657.6928567215408
Iteration:      5,   Func. Count:     46,   Neg. LLF: -657.7021169668074
Iteration:      6,   Func. Count:     55,   Neg. LLF: -657.7205035925924
Iteration:      7,   Func. Count:     63,   Neg. LLF: -657.7315250892177
Iteration:      8,   Func. Count:     72,   Neg. LLF: -657.7330213361397
Iteration:      9,   Func. Count:     79,   Neg. LLF: -657.7341256243355
Iteration:     10,   Func. Count:     86,   Neg. LLF: -657.7342081525653
Iteration:     11,   Func. Count:     91,   Neg. LLF: -657.7342209027067
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -657.7342209027092
            Iterations: 12
            Function evaluations: 91
            Gradient evaluations: 11
Iteration:      1,   Func. Count:      5,   Neg. LLF: -655.9892348458227
Iteration:      2,   Func. Count:     17,   Neg. LLF: -656.0590018942496
Iteration:      3,   Func. Count:     26,   Neg. LLF: -656.0642219327943
Iteration:      4,   Func. Count:     38,   Neg. LLF: -656.0681600579908
Iteration:      5,   Func. Count:     46,   Neg. LLF: -656.0774570732628
Iteration:      6,   Func. Count:     55,   Neg. LLF: -656.0790705839521
Iteration:      7,   Func. Count:     62,   Neg. LLF: -656.0979807504824
Iteration:      8,   Func. Count:     71,   Neg. LLF: -656.0979892332143
Iteration:      9,   Func. Count:     76,   Neg. LLF: -656.0980105987001
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -656.0980105987187
            Iterations: 10
            Function evaluations: 76
            Gradient evaluations: 9
Iteration:      1,   Func. Count:      5,   Neg. LLF: -655.9042540899916
Iteration:      2,   Func. Count:     17,   Neg. LLF: -655.9721247197159
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -655.9721247101422
            Iterations: 3
            Function evaluations: 26
            Gradient evaluations: 2
Iteration:      1,   Func. Count:      5,   Neg. LLF: -654.4405482551299
Iteration:      2,   Func. Count:     17,   Neg. LLF: -654.4851644066307
Iteration:      3,   Func. Count:     25,   Neg. LLF: -654.4894933490531
Iteration:      4,   Func. Count:     32,   Neg. LLF: -654.5311277541294
Iteration:      5,   Func. Count:     41,   Neg. LLF: -654.5312516979479
Iteration:      6,   Func. Count:     46,   Neg. LLF: -654.5312592957482
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -654.5312592957519
            Iterations: 6
            Function evaluations: 46
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -654.5182585181257
Iteration:      2,   Func. Count:     17,   Neg. LLF: -654.5657937069122
Iteration:      3,   Func. Count:     25,   Neg. LLF: -654.5663402450866
Iteration:      4,   Func. Count:     37,   Neg. LLF: -654.5668412610971
Iteration:      5,   Func. Count:     44,   Neg. LLF: -654.5727006585873
Iteration:      6,   Func. Count:     55,   Neg. LLF: -654.5727028954988
Iteration:      7,   Func. Count:     67,   Neg. LLF: -654.5727077336569
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -654.5727077598408
            Iterations: 9
            Function evaluations: 74
            Gradient evaluations: 7
Iteration:      1,   Func. Count:      5,   Neg. LLF: -651.1112329638283
Iteration:      2,   Func. Count:     17,   Neg. LLF: -651.1331228514488
Iteration:      3,   Func. Count:     25,   Neg. LLF: -651.1341130783869
Iteration:      4,   Func. Count:     32,   Neg. LLF: -651.152273129842
Iteration:      5,   Func. Count:     40,   Neg. LLF: -651.152648565452
Iteration:      6,   Func. Count:     48,   Neg. LLF: -651.153704766913
Iteration:      7,   Func. Count:     53,   Neg. LLF: -651.1540469025593
Iteration:      8,   Func. Count:     58,   Neg. LLF: -651.1540645421729
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -651.1540645421821
            Iterations: 8
            Function evaluations: 58
            Gradient evaluations: 8
Iteration:      1,   Func. Count:      5,   Neg. LLF: -649.204720775396
Iteration:      2,   Func. Count:     17,   Neg. LLF: -649.2698276097943
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -649.2698280637319
            Iterations: 2
            Function evaluations: 23
            Gradient evaluations: 2
Iteration:      1,   Func. Count:      5,   Neg. LLF: -647.8648840508957
Iteration:      2,   Func. Count:     17,   Neg. LLF: -647.9344547181375
Iteration:      3,   Func. Count:     25,   Neg. LLF: -647.9419209023314
Iteration:      4,   Func. Count:     32,   Neg. LLF: -647.9497555301346
Iteration:      5,   Func. Count:     39,   Neg. LLF: -647.9648216657733
Iteration:      6,   Func. Count:     47,   Neg. LLF: -647.9648476529103
Iteration:      7,   Func. Count:     52,   Neg. LLF: -647.9648619540646
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -647.9648619540735
            Iterations: 7
            Function evaluations: 52
            Gradient evaluations: 7
Iteration:      1,   Func. Count:      5,   Neg. LLF: -651.4129660627998
Iteration:      2,   Func. Count:     17,   Neg. LLF: -651.4575889595142
Iteration:      3,   Func. Count:     28,   Neg. LLF: -651.4575904571143
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -651.4575904335206
            Iterations: 5
            Function evaluations: 37
            Gradient evaluations: 3
Iteration:      1,   Func. Count:      5,   Neg. LLF: -649.3859898771539
Iteration:      2,   Func. Count:     17,   Neg. LLF: -649.5104122355534
Iteration:      3,   Func. Count:     26,   Neg. LLF: -649.5108472356878
Iteration:      4,   Func. Count:     33,   Neg. LLF: -649.5187741558833
Iteration:      5,   Func. Count:     41,   Neg. LLF: -649.5187779292498
Iteration:      6,   Func. Count:     46,   Neg. LLF: -649.5188578227564
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -649.5188578227702
            Iterations: 6
            Function evaluations: 46
            Gradient evaluations: 6
Iteration:      1,   Func. Count:      5,   Neg. LLF: -649.4320671307988
Iteration:      2,   Func. Count:     17,   Neg. LLF: -649.5770448825115
Iteration:      3,   Func. Count:     26,   Neg. LLF: -649.577118372875
Iteration:      4,   Func. Count:     33,   Neg. LLF: -649.578102911356
Iteration:      5,   Func. Count:     38,   Neg. LLF: -649.5781516380657
Iteration:      6,   Func. Count:     44,   Neg. LLF: -649.578230463566
Iteration:      7,   Func. Count:     49,   Neg. LLF: -649.5782586700102
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -649.5782586700732
            Iterations: 7
            Function evaluations: 49
            Gradient evaluations: 7
Iteration:      1,   Func. Count:      5,   Neg. LLF: -649.8110083626044
Iteration:      2,   Func. Count:     17,   Neg. LLF: -649.9318119205432
Iteration:      3,   Func. Count:     29,   Neg. LLF: -649.9318146518306
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -649.9318150977472
            Iterations: 5
            Function evaluations: 37
            Gradient evaluations: 3
Iteration:      1,   Func. Count:      5,   Neg. LLF: -650.805578132557
Iteration:      2,   Func. Count:     17,   Neg. LLF: -650.9326486836682
Iteration:      3,   Func. Count:     25,   Neg. LLF: -650.937820495724
Iteration:      4,   Func. Count:     33,   Neg. LLF: -650.9384003660639
Iteration:      5,   Func. Count:     39,   Neg. LLF: -650.9419266747137
Iteration:      6,   Func. Count:     48,   Neg. LLF: -650.9419372375771
Iteration:      7,   Func. Count:     54,   Neg. LLF: -650.946890739448
Iteration:      8,   Func. Count:     59,   Neg. LLF: -650.9468919094435
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -650.9468919094534
            Iterations: 8
            Function evaluations: 59
            Gradient evaluations: 8
Iteration:      1,   Func. Count:      5,   Neg. LLF: -650.4885855609117
Iteration:      2,   Func. Count:     17,   Neg. LLF: -650.5533213306016
Iteration:      3,   Func. Count:     26,   Neg. LLF: -650.553568046112
Iteration:      4,   Func. Count:     33,   Neg. LLF: -650.5539672212103
Iteration:      5,   Func. Count:     38,   Neg. LLF: -650.554126079099
Iteration:      6,   Func. Count:     44,   Neg. LLF: -650.5543007638967
Iteration:      7,   Func. Count:     49,   Neg. LLF: -650.5543643864472
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -650.5543643865449
            Iterations: 7
            Function evaluations: 49
            Gradient evaluations: 7
Iteration:      1,   Func. Count:      5,   Neg. LLF: -655.3122475582609
Iteration:      2,   Func. Count:     17,   Neg. LLF: -655.3419618895005
Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: -655.3419618740463
            Iterations: 6
            Function evaluations: 17
            Gradient evaluations: 2

/Users/austincostello1/anaconda3/lib/python3.6/site-packages/arch/univariate/base.py:522: ConvergenceWarning: 
The optimizer returned code 8. The message is:
Positive directional derivative for linesearch
See scipy.optimize.fmin_slsqp for code meaning.

  ConvergenceWarning)
/Users/austincostello1/anaconda3/lib/python3.6/site-packages/arch/univariate/base.py:522: ConvergenceWarning: 
The optimizer returned code 8. The message is:
Positive directional derivative for linesearch
See scipy.optimize.fmin_slsqp for code meaning.

  ConvergenceWarning)
/Users/austincostello1/anaconda3/lib/python3.6/site-packages/arch/univariate/base.py:522: ConvergenceWarning: 
The optimizer returned code 8. The message is:
Positive directional derivative for linesearch
See scipy.optimize.fmin_slsqp for code meaning.

  ConvergenceWarning)

PRED_1 = np.reshape(PRED, (589,6))

Let's check how the model is performing:

PRED_1[:,0]
plt.plot(ACTUAL, label='Actual')
plt.plot(PRED_1[:,0], label='Predicted')
plt.title("SPY Variance vs Predicted Variance from ARCH(1)")
plt.legend()
plt.savefig("arch.png", bbox_inches='tight')

rolling_ar_errors = measure_error(ACTUAL, PRED_1[:,0], label='RMSE ARCH(1)'  )

RMSE = 0.00015760018121505748
RMSE_pcent = 1.9817494383790915
label = RMSE ARCH(1)

Time to check if we made money with this strategy:¶

## Calculating PNL 
dir = []
for i in range(len(PRED_1)-1):
    for j in range(len(PRED_1[0])):
        if ACTUAL[i] < PRED_1[i+1][j]:
            dir.append(1)
        else:
            dir.append(-1)

dir1=np.reshape(dir, (588, 6))

act_vix_lst = df_data.loc[df_data.Date > '2014-07-29',["VIX","VIXCM30","VIXCM60","VIXCM90","VIXCM120","VIXCM150"]]
act_vix_lst = act_vix_lst.values
diff = np.diff(act_vix_lst, axis=0)

money_accumulated = np.multiply(diff,dir1)
PnL = np.sum(money_accumulated)
print(PnL)

-58.35501325999999

Summary of ARCH and GARCH Models¶

model	PnL	RMSE
ARCH(1)	-58.3550	0.000157
GARCH(1,1)	-78.4006	0.000153
GARCH(2,1)	-99.9630	0.000152
GARCH(2,2)	-91.0460	0.000150
GARCH(3,2)	-95.5904	0.000152

While our RMSEs look quite strong, none of these models seem to result in a profitable strategy, so we will move onto try VECM models.

VECM
¶

VIXCM30
¶

VIXCM30¶

First let's take a look at the plotted data:

df_data.VIX.plot(label='VIX')
df_data.VIXCM30.plot(label='VIXCM30')
plt.legend()

<matplotlib.legend.Legend at 0x1c1dc3e160>

The two variables look cointegrated, but let's run a Johansen Test to be sure:

d = {'vix': df_data['VIX'],
     'vixcm30':df_data['VIXCM30']
    }
df = pd.DataFrame(d)
c = coint_johansen(df,1,63)
c

--------------------------------------------------
--> Trace Statistics
variable statistic Crit-90% Crit-95%  Crit-99%
r = 0 	 39.5184 16.1619 18.3985 23.1485
r = 1 	 8.231 2.7055 3.8415 6.6349
--------------------------------------------------
--> Eigen Statistics
variable statistic Crit-90% Crit-95%  Crit-99%
r = 0 	 31.2874 15.0006 17.1481 21.7465
r = 1 	 8.231 2.7055 3.8415 6.6349
--------------------------------------------------
eigenvectors:
 [[ 0.77065652 -0.08729355]
 [-0.83501063  0.23909376]]
--------------------------------------------------
eigenvalues:
 [0.0108086 0.0028549]
--------------------------------------------------

<johansen_test.Holder at 0x1c1d273940>

The test confirms they are cointegrated so we can build the model:

train,test = train_test_split(df_data, test_size=0.2, shuffle=False)

H = train[['VIX', 'VIXCM30']].dropna()
v = VECM(H, exog=train['SPY'],deterministic='li', k_ar_diff=20)
f = v.fit()
f.summary()

Picking up residuals for VIXCM30 series from the model:

# getting the lag to remove residual information
resid_vixcm30 = []
for i in range(len(f.resid)):
    resid_vixcm30 += [f.resid[i][1]]

plt.plot(resid_vixcm30)

[<matplotlib.lines.Line2D at 0x1c2392dc50>]

run_ljung_box(resid_vixcm30, 80)

Ljung-Box Test

No information is left in the residuals!

Now we can check how the model is performing:

H = train[['VIX', 'VIXCM30', 'SPY']].dropna()
res = []
for date, row in test.iterrows(): 
    v = VECM(H[['VIX', 'VIXCM30']], exog=H['SPY'], deterministic='li', k_ar_diff=20)
    f = v.fit()
    res += [f.predict(steps=1, exog_fc=test['SPY'])[0][1]]
    H = H.append(row[['VIX', 'VIXCM30', 'SPY']])
    H = H[1:].dropna()

plot_df = pd.DataFrame()
plot_df['Date'] = test.iloc[:, 0].values
plot_df['VIXCM30_actual'] = test.iloc[:, 4].values
plot_df['VIXCM30_predicted'] = res

plot_df.VIXCM30_actual.plot()
plot_df.VIXCM30_predicted.plot()
plt.title("VIXCM30 prediction through VECM model")
plt.legend()

<matplotlib.legend.Legend at 0x1c23938b38>

Looking closely into a random sample of the predictions:

plot_df.VIXCM30_actual[-100:].plot()
plot_df.VIXCM30_predicted[-100:].plot()
plt.title("VIXCM30 prediction through VECM model")
plt.legend()

<matplotlib.legend.Legend at 0x1c24364ba8>

ar_error2 = measure_error(plot_df['VIXCM30_actual'].values, plot_df['VIXCM30_predicted'].values, label='VECM_30')
print(ar_error2)

RMSE = 0.8218564880183119
RMSE_pcent = 0.04658378560777323
label = VECM_30
{'RMSE': 0.8218564880183119, 'RMSE_pcent': 0.04658378560777323, 'label': 'VECM_30'}

make_money(plot_df.VIXCM30_actual, plot_df.VIXCM30_predicted)

Total profit or loss: 17.6314
Number of correct days: 310
total incorrect days: 278

VIXCM60
¶

VIXCM60¶

First let's take a look at the plotted data:

df_data.VIX.plot(label='VIX')
df_data.VIXCM60.plot(label='VIXCM60')
plt.legend()

<matplotlib.legend.Legend at 0x1c20485e48>

The two variables look cointegrated, but let's run a Johansen Test to be sure:

d = {'vix': df_data['VIX'],
     'vixcm60':df_data['VIXCM60']
    }
df = pd.DataFrame(d)
c = coint_johansen(df,1,63)
c

--------------------------------------------------
--> Trace Statistics
variable statistic Crit-90% Crit-95%  Crit-99%
r = 0 	 30.388 16.1619 18.3985 23.1485
r = 1 	 7.2609 2.7055 3.8415 6.6349
--------------------------------------------------
--> Eigen Statistics
variable statistic Crit-90% Crit-95%  Crit-99%
r = 0 	 23.1271 15.0006 17.1481 21.7465
r = 1 	 7.2609 2.7055 3.8415 6.6349
--------------------------------------------------
eigenvectors:
 [[ 0.48077253 -0.10256322]
 [-0.52658657  0.27381262]]
--------------------------------------------------
eigenvalues:
 [0.00800084 0.00251885]
--------------------------------------------------

<johansen_test.Holder at 0x113008c50>

Now we can build the model:

H = train[['VIX', 'VIXCM60']].dropna()
v = VECM(H, exog=train['SPY'],deterministic='li', k_ar_diff=20)
f = v.fit()
f.summary()

# getting the lag to remove residual information
resid_vixcm60 = []
for i in range(len(f.resid)):
    resid_vixcm60 += [f.resid[i][1]]

plt.plot(resid_vixcm60)

[<matplotlib.lines.Line2D at 0x1c269f4c88>]

run_ljung_box(resid_vixcm60, 80)

Ljung-Box Test

Now we can check how the model is performing:

H = train[['VIX', 'VIXCM60', 'SPY']].dropna()
res = []
for date, row in test.iterrows(): 
    v = VECM(H[['VIX', 'VIXCM60']], exog=H['SPY'], deterministic='li', k_ar_diff=20)
    f = v.fit()
    res += [f.predict(steps=1, exog_fc=test['SPY'])[0][1]]
    H = H.append(row[['VIX', 'VIXCM60', 'SPY']])
    H = H[1:].dropna()

plot_df = pd.DataFrame()
plot_df['Date'] = test.iloc[:, 0].values
plot_df['VIXCM60_actual'] = test.iloc[:, 5].values
plot_df['VIXCM60_predicted'] = res

plot_df.VIXCM60_actual.plot()
plot_df.VIXCM60_predicted.plot()
plt.title("VIXCM60 prediction through VECM model")
plt.legend()

<matplotlib.legend.Legend at 0x1c26cf44a8>

ar_error2 = measure_error(plot_df['VIXCM60_actual'].values, plot_df['VIXCM60_predicted'].values, label='VECM_60')
print(ar_error2)

RMSE = 0.619971615010414
RMSE_pcent = 0.03360636343771486
label = VECM_60
{'RMSE': 0.619971615010414, 'RMSE_pcent': 0.03360636343771486, 'label': 'VECM_60'}

make_money(plot_df.VIXCM60_actual, plot_df.VIXCM60_predicted)

Total profit or loss: -8.6022
Number of correct days: 301
total incorrect days: 287

VIXCM90
¶

VIXCM90¶

First let's take a look at the plotted data:

df_data.VIX.plot(label='VIX')
df_data.VIXCM90.plot(label='VIXCM90')
plt.legend()

<matplotlib.legend.Legend at 0x1c26d31908>

The two variables look cointegrated, but let's run a Johansen Test to be sure:

d = {'vix': df_data['VIX'],
     'vixcm90':df_data['VIXCM90']
    }
df = pd.DataFrame(d)
c = coint_johansen(df,1,63)
c

--------------------------------------------------
--> Trace Statistics
variable statistic Crit-90% Crit-95%  Crit-99%
r = 0 	 30.0345 16.1619 18.3985 23.1485
r = 1 	 7.0617 2.7055 3.8415 6.6349
--------------------------------------------------
--> Eigen Statistics
variable statistic Crit-90% Crit-95%  Crit-99%
r = 0 	 22.9729 15.0006 17.1481 21.7465
r = 1 	 7.0617 2.7055 3.8415 6.6349
--------------------------------------------------
eigenvectors:
 [[ 0.38643582 -0.07550673]
 [-0.42911003  0.25519026]]
--------------------------------------------------
eigenvalues:
 [0.00794771 0.00244981]
--------------------------------------------------

<johansen_test.Holder at 0x1c26db4dd8>

Now we can build the model:

H = train[['VIX', 'VIXCM90']].dropna()
v = VECM(H, deterministic='li',exog=train['SPY'], k_ar_diff=20)
f = v.fit()
f.summary()

# getting the lag to remove residual information
resid_vixcm90 = []
for i in range(len(f.resid)):
    resid_vixcm90 += [f.resid[i][1]]

plt.plot(resid_vixcm90)

[<matplotlib.lines.Line2D at 0x1c2772deb8>]

run_ljung_box(resid_vixcm90, 80)

Ljung-Box Test

Now we can check how the model is performing:

H = train[['VIX', 'VIXCM90', 'SPY']].dropna()
res = []
for date, row in test.iterrows(): 
    v = VECM(H[['VIX', 'VIXCM90']], exog=H['SPY'], deterministic='li', k_ar_diff=20)
    f = v.fit()
    res += [f.predict(steps=1, exog_fc=test['SPY'])[0][1]]
    H = H.append(row[['VIX', 'VIXCM90', 'SPY']])
    H = H[1:].dropna()

plot_df = pd.DataFrame()
plot_df['Date'] = test.iloc[:,0].values
plot_df['VIXCM90_actual'] = test.iloc[:, 6].values
plot_df['VIXCM90_predicted'] = res

plot_df.VIXCM90_actual.plot()
plot_df.VIXCM90_predicted.plot()
plt.title("VIXCM90 prediction through VECM model")
plt.legend()

<matplotlib.legend.Legend at 0x1c278c5400>

ar_error2 = measure_error(plot_df['VIXCM90_actual'].values, plot_df['VIXCM90_predicted'].values, label='VECM_90')
print(ar_error2)

RMSE = 0.5184302696563623
RMSE_pcent = 0.027288800658595697
label = VECM_90
{'RMSE': 0.5184302696563623, 'RMSE_pcent': 0.027288800658595697, 'label': 'VECM_90'}

make_money(plot_df.VIXCM90_actual, plot_df.VIXCM90_predicted)

Total profit or loss: -13.8828
Number of correct days: 286
total incorrect days: 302

VIXCM120
¶

VIXCM120¶

First let's take a look at the plotted data:

df_data.VIX.plot(label='VIX')
df_data.VIXCM120.plot(label='VIXCM120')
plt.legend()

<matplotlib.legend.Legend at 0x1c27e4de10>

The two variables look cointegrated, but let's run a Johansen Test to be sure:

d = {'vix': df_data['VIX'],
     'vixcm120':df_data['VIXCM120']
    }
df = pd.DataFrame(d)
c = coint_johansen(df,1,63)
c

--------------------------------------------------
--> Trace Statistics
variable statistic Crit-90% Crit-95%  Crit-99%
r = 0 	 26.3306 16.1619 18.3985 23.1485
r = 1 	 6.8779 2.7055 3.8415 6.6349
--------------------------------------------------
--> Eigen Statistics
variable statistic Crit-90% Crit-95%  Crit-99%
r = 0 	 19.4527 15.0006 17.1481 21.7465
r = 1 	 6.8779 2.7055 3.8415 6.6349
--------------------------------------------------
eigenvectors:
 [[ 0.3257479  -0.08536839]
 [-0.35188173  0.27314798]]
--------------------------------------------------
eigenvalues:
 [0.00673399 0.00238613]
--------------------------------------------------

<johansen_test.Holder at 0x10bdc7e10>

Now we can build the model:

H = train[['VIX', 'VIXCM120']].dropna()
v = VECM(H, deterministic='li', exog=train['SPY'], k_ar_diff=20)
f = v.fit()
f.summary()

# getting the lag to remove residual information
resid_vixcm120 = []
for i in range(len(f.resid)):
    resid_vixcm120 += [f.resid[i][1]]

plt.plot(resid_vixcm120)

[<matplotlib.lines.Line2D at 0x1c204d5cc0>]

run_ljung_box(resid_vixcm120, 80)

Ljung-Box Test

Now we can check how the model is performing:

H = train[['VIX', 'VIXCM120', 'SPY']].dropna()
res = []
for date, row in test.iterrows(): 
    v = VECM(H[['VIX', 'VIXCM120']], exog=H['SPY'], deterministic='li', k_ar_diff=20)
    f = v.fit()
    res += [f.predict(steps=1, exog_fc=test['SPY'])[0][1]]
    H = H.append(row[['VIX', 'VIXCM120', 'SPY']])
    H = H[1:].dropna()

plot_df = pd.DataFrame()
plot_df['Date'] = test.iloc[:,0].values
plot_df['VIXCM120_actual'] = test.iloc[:,7].values
plot_df['VIXCM120_predicted'] = res

plot_df.VIXCM120_actual.plot()
plot_df.VIXCM120_predicted.plot()
plt.title("VIXCM120 prediction through VECM model")
plt.legend()

<matplotlib.legend.Legend at 0x1c27983eb8>

ar_error2 = measure_error(plot_df['VIXCM120_actual'].values, plot_df['VIXCM120_predicted'].values, label='VECM')
print(ar_error2)

RMSE = 0.4581289208624565
RMSE_pcent = 0.023545101076493956
label = VECM
{'RMSE': 0.4581289208624565, 'RMSE_pcent': 0.023545101076493956, 'label': 'VECM'}

make_money(plot_df.VIXCM120_actual, plot_df.VIXCM120_predicted)

Total profit or loss: -7.5644
Number of correct days: 294
total incorrect days: 294

VIXCM150
¶

VIXCM150¶

First let's take a look at the plotted data:

df_data.VIX.plot(label='VIX')
df_data.VIXCM150.plot(label='VIXCM150')
plt.legend()

<matplotlib.legend.Legend at 0x1c204a8fd0>

The two variables look cointegrated, but let's run a Johansen Test to be sure:

d = {'vix': df_data['VIX'],
     'vixcm150':df_data['VIXCM150']
    }
df = pd.DataFrame(d)
c = coint_johansen(df,1,63)
c

--------------------------------------------------
--> Trace Statistics
variable statistic Crit-90% Crit-95%  Crit-99%
r = 0 	 25.2241 16.1619 18.3985 23.1485
r = 1 	 6.6399 2.7055 3.8415 6.6349
--------------------------------------------------
--> Eigen Statistics
variable statistic Crit-90% Crit-95%  Crit-99%
r = 0 	 18.5842 15.0006 17.1481 21.7465
r = 1 	 6.6399 2.7055 3.8415 6.6349
--------------------------------------------------
eigenvectors:
 [[ 0.2736747  -0.08686034]
 [-0.28087418  0.27910247]]
--------------------------------------------------
eigenvalues:
 [0.00643429 0.00230366]
--------------------------------------------------

<johansen_test.Holder at 0x1c27c45470>

Now we can build the model:

H = train[['VIX', 'VIXCM150']].dropna()
v = VECM(H, deterministic='li', exog=train['SPY'], k_ar_diff=20)
f = v.fit()
f.summary()

# getting the lag to remove residual information
resid_vixcm150 = []
for i in range(len(f.resid)):
    resid_vixcm150 += [f.resid[i][1]]

plt.plot(resid_vixcm150)

[<matplotlib.lines.Line2D at 0x1c1ccde358>]

run_ljung_box(resid_vixcm150, 80)

Ljung-Box Test

Now we can check how the model is performing:

H = train[['VIX', 'VIXCM150', 'SPY']].dropna()
res = []
for date, row in test.iterrows(): 
    v = VECM(H[['VIX', 'VIXCM150']], exog=H['SPY'], deterministic='li', k_ar_diff=20)
    f = v.fit()
    res += [f.predict(steps=1, exog_fc=test['SPY'])[0][1]]
    H = H.append(row[['VIX', 'VIXCM150', 'SPY']])
    H = H[1:].dropna()

plot_df = pd.DataFrame()
plot_df['Date'] = test.iloc[:, 0].values
plot_df['VIXCM150_actual'] = test.iloc[:,8].values
plot_df['VIXCM150_predicted'] = res

plot_df.VIXCM150_actual.plot()
plot_df.VIXCM150_predicted.plot()
plt.title("VIXCM150 prediction through VECM model")
plt.legend()

<matplotlib.legend.Legend at 0x1c27cc8748>

ar_error2 = measure_error(plot_df['VIXCM150_actual'].values, plot_df['VIXCM150_predicted'].values, label='VECM')
print(ar_error2)

RMSE = 0.39943060263751484
RMSE_pcent = 0.020057050701287188
label = VECM
{'RMSE': 0.39943060263751484, 'RMSE_pcent': 0.020057050701287188, 'label': 'VECM'}

make_money(plot_df.VIXCM150_actual, plot_df.VIXCM150_predicted)

Total profit or loss: -5.7534
Number of correct days: 290
total incorrect days: 298

Summary of VECM Models
¶

Summary of VECM Models:¶

Model	RMSE	PnL
VIXCM30	0.8218564	17.6314
VIXCM60	0.6199716	-8.6022
VIXCM90	0.5184302	-13.8828
VIXCM120	0.4581289	-7.5644
VIXCM150	0.3994306	-5.7534

AR Models
¶

VIXCM30 Premium
VIXCM60 Premium
VIXCM90 Premium
VIXCM120 Premium
VIXCM150 Premium
Summary of AR Models

Defining New Premium Variables:¶

For our AR and ARMA models we will try to forecast the premium of each VIX future. Since the price of these futures is the VIX + a premium (cost of carry) we can create a variable for the premium by subtracting the VIX from the price of any given future. We do this for each future below. After running the model on these premium prices, we will add the most recent VIX price back to the premium to find the VIXCM future price for the next day.

df_data.set_index('Date', inplace= True)

df_data.head()

df_data['VIXCM30_Premium'] = (df_data['VIXCM30']-df_data['VIX'])

df_data['VIXCM60_Premium'] = (df_data['VIXCM60']-df_data['VIX'])

df_data['VIXCM90_Premium'] = (df_data['VIXCM90']-df_data['VIX'])

df_data['VIXCM120_Premium'] = (df_data['VIXCM120']-df_data['VIX'])

df_data['VIXCM150_Premium'] = (df_data['VIXCM150']-df_data['VIX'])

Now we will check if all of the premium variables are stationary:¶

VIXCM30 Premium

df_data.VIXCM30_Premium.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c19c1a828>

run_ad_fuller(df_data.VIXCM30_Premium)

ADF Statistic: -6.546746
p-value: 0.000000
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567

We reject the Null Hypothesis (H0) -- the time series is stationary

It is Stationary!

VIXCM60 Premium

df_data.VIXCM60_Premium.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c19f4fa20>

run_ad_fuller(df_data.VIXCM60_Premium)

ADF Statistic: -5.490771
p-value: 0.000002
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567

We reject the Null Hypothesis (H0) -- the time series is stationary

It is Stationary!

VIXCM90 Premium

df_data.VIXCM90_Premium.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x10acb0a20>

run_ad_fuller(df_data.VIXCM90_Premium)

ADF Statistic: -5.085105
p-value: 0.000015
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567

We reject the Null Hypothesis (H0) -- the time series is stationary

It is Stationary!

VIXCM120 Premium

df_data.VIXCM120_Premium.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c19d34cc0>

run_ad_fuller(df_data.VIXCM120_Premium.dropna())

ADF Statistic: -4.742800
p-value: 0.000070
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567

We reject the Null Hypothesis (H0) -- the time series is stationary

It is Stationary!

VIXCM150 Premium

df_data.VIXCM150_Premium.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c1d4eb908>

run_ad_fuller(df_data.VIXCM150_Premium.dropna())

ADF Statistic: -4.555457
p-value: 0.000156
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567

We reject the Null Hypothesis (H0) -- the time series is stationary

It is Stationary!

VIXCM30 Premium
¶

AR Model for VIXCM30 Premium¶

First we will split the data into test and training data:

train,test = train_test_split(df_data, test_size=0.2, shuffle=False)
print('train: ',len(train))
print('test : ',len(test))
train.tail(5)

train:  2354
test :  589

Now we will look at the PACF plot of the premimums:

plot_pacf(train['VIXCM30_Premium'],lags=50)
plt.show()

# train autoregression
ar_model = AR(train['VIXCM30_Premium'])

Now we will check what order model is best from AIC and BIC scores:

aic = ar_model.select_order(maxlag=30,ic='aic',trend='nc')
bic = ar_model.select_order(maxlag=30,ic='bic',trend='nc')

print('model order (AIC): %d' % aic)
print('model order (BIC): %d' % bic)

model order (AIC): 30
model order (BIC): 8

Above we can see that the AIC score suggests 30 lags is the best, but this is the max lag we allowed it to use, so we will increase the max to see if it tells us something different.

aic = ar_model.select_order(maxlag=50,ic='aic',trend='nc')
bic = ar_model.select_order(maxlag=50,ic='bic',trend='nc')

print('model order (AIC): %d' % aic)
print('model order (BIC): %d' % bic)

model order (AIC): 49
model order (BIC): 8

We will try order =8 since it is the simplest model:

model_fit = ar_model.fit(maxlag=8)
print('Lag: %s' % model_fit.k_ar)
print('Coefficients: %s' % model_fit.params)
intercept = model_fit.params[0]
coeff = list(model_fit.params)[1:] # drop the intercept
coeff.reverse()

Lag: 8
Coefficients: const                 0.040676
L1.VIXCM30_Premium    0.739705
L2.VIXCM30_Premium    0.066149
L3.VIXCM30_Premium    0.052504
L4.VIXCM30_Premium   -0.003647
L5.VIXCM30_Premium    0.110634
L6.VIXCM30_Premium   -0.049825
L7.VIXCM30_Premium   -0.042842
L8.VIXCM30_Premium    0.089621
dtype: float64

print('RMSE (in sample):',np.sqrt(np.mean(model_fit.sigma2)))

RMSE (in sample): 1.191782092547236

We will check to make sure that there is no information left in the residuals and that they are stationary:

run_ljung_box(model_fit.resid, lags=100)

Ljung-Box Test
lag:  10 | p_value: 0.0275 | still has information
lag:  11 | p_value: 0.0354 | still has information
lag:  12 | p_value: 0.0318 | still has information
lag:  13 | p_value: 0.0450 | still has information
lag:  14 | p_value: 0.0258 | still has information
lag:  15 | p_value: 0.0349 | still has information
lag:  16 | p_value: 0.0264 | still has information
lag:  17 | p_value: 0.0330 | still has information
lag:  18 | p_value: 0.0098 | still has information
lag:  19 | p_value: 0.0072 | still has information
lag:  20 | p_value: 0.0106 | still has information
lag:  21 | p_value: 0.0058 | still has information
lag:  22 | p_value: 0.0078 | still has information
lag:  23 | p_value: 0.0072 | still has information
lag:  24 | p_value: 0.0072 | still has information
lag:  25 | p_value: 0.0090 | still has information
lag:  26 | p_value: 0.0087 | still has information
lag:  27 | p_value: 0.0035 | still has information
lag:  28 | p_value: 0.0049 | still has information
lag:  29 | p_value: 0.0008 | still has information
lag:  30 | p_value: 0.0005 | still has information
lag:  31 | p_value: 0.0004 | still has information
lag:  32 | p_value: 0.0005 | still has information
lag:  33 | p_value: 0.0001 | still has information
lag:  34 | p_value: 0.0000 | still has information
lag:  35 | p_value: 0.0000 | still has information
lag:  36 | p_value: 0.0000 | still has information
lag:  37 | p_value: 0.0000 | still has information
lag:  38 | p_value: 0.0000 | still has information
lag:  39 | p_value: 0.0000 | still has information
lag:  40 | p_value: 0.0000 | still has information
lag:  41 | p_value: 0.0000 | still has information
lag:  42 | p_value: 0.0000 | still has information
lag:  43 | p_value: 0.0000 | still has information
lag:  44 | p_value: 0.0000 | still has information
lag:  45 | p_value: 0.0000 | still has information
lag:  46 | p_value: 0.0000 | still has information
lag:  47 | p_value: 0.0000 | still has information
lag:  48 | p_value: 0.0000 | still has information
lag:  49 | p_value: 0.0000 | still has information
lag:  50 | p_value: 0.0000 | still has information
lag:  51 | p_value: 0.0000 | still has information
lag:  52 | p_value: 0.0000 | still has information
lag:  53 | p_value: 0.0000 | still has information
lag:  54 | p_value: 0.0000 | still has information
lag:  55 | p_value: 0.0000 | still has information
lag:  56 | p_value: 0.0000 | still has information
lag:  57 | p_value: 0.0000 | still has information
lag:  58 | p_value: 0.0000 | still has information
lag:  59 | p_value: 0.0000 | still has information
lag:  60 | p_value: 0.0000 | still has information
lag:  61 | p_value: 0.0000 | still has information
lag:  62 | p_value: 0.0000 | still has information
lag:  63 | p_value: 0.0000 | still has information
lag:  64 | p_value: 0.0000 | still has information
lag:  65 | p_value: 0.0000 | still has information
lag:  66 | p_value: 0.0000 | still has information
lag:  67 | p_value: 0.0000 | still has information
lag:  68 | p_value: 0.0000 | still has information
lag:  69 | p_value: 0.0000 | still has information
lag:  70 | p_value: 0.0000 | still has information
lag:  71 | p_value: 0.0000 | still has information
lag:  72 | p_value: 0.0000 | still has information
lag:  73 | p_value: 0.0000 | still has information
lag:  74 | p_value: 0.0000 | still has information
lag:  75 | p_value: 0.0000 | still has information
lag:  76 | p_value: 0.0000 | still has information
lag:  77 | p_value: 0.0000 | still has information
lag:  78 | p_value: 0.0000 | still has information
lag:  79 | p_value: 0.0000 | still has information
lag:  80 | p_value: 0.0000 | still has information
lag:  81 | p_value: 0.0000 | still has information
lag:  82 | p_value: 0.0000 | still has information
lag:  83 | p_value: 0.0000 | still has information
lag:  84 | p_value: 0.0000 | still has information
lag:  85 | p_value: 0.0000 | still has information
lag:  86 | p_value: 0.0000 | still has information
lag:  87 | p_value: 0.0000 | still has information
lag:  88 | p_value: 0.0000 | still has information
lag:  89 | p_value: 0.0000 | still has information
lag:  90 | p_value: 0.0000 | still has information
lag:  91 | p_value: 0.0000 | still has information
lag:  92 | p_value: 0.0000 | still has information
lag:  93 | p_value: 0.0000 | still has information
lag:  94 | p_value: 0.0000 | still has information
lag:  95 | p_value: 0.0000 | still has information
lag:  96 | p_value: 0.0000 | still has information
lag:  97 | p_value: 0.0000 | still has information
lag:  98 | p_value: 0.0000 | still has information
lag:  99 | p_value: 0.0000 | still has information
lag: 100 | p_value: 0.0000 | still has information

Now we will try order =49 to see if we can caputure this missing information:

model_fit = ar_model.fit(maxlag=49)
print('Lag: %s' % model_fit.k_ar)
print('Coefficients: %s' % model_fit.params)
intercept = model_fit.params[0]
coeff = list(model_fit.params)[1:] # drop the intercept
coeff.reverse()

Lag: 49
Coefficients: const                  0.047729
L1.VIXCM30_Premium     0.734760
L2.VIXCM30_Premium     0.073003
L3.VIXCM30_Premium     0.054725
L4.VIXCM30_Premium     0.006517
L5.VIXCM30_Premium     0.089636
L6.VIXCM30_Premium    -0.042873
L7.VIXCM30_Premium    -0.024546
L8.VIXCM30_Premium     0.086363
L9.VIXCM30_Premium     0.006226
L10.VIXCM30_Premium    0.077364
L11.VIXCM30_Premium   -0.058672
L12.VIXCM30_Premium   -0.046534
L13.VIXCM30_Premium    0.010955
L14.VIXCM30_Premium    0.045101
L15.VIXCM30_Premium   -0.023944
L16.VIXCM30_Premium    0.019762
L17.VIXCM30_Premium   -0.010226
L18.VIXCM30_Premium   -0.065961
L19.VIXCM30_Premium    0.061991
L20.VIXCM30_Premium   -0.034029
L21.VIXCM30_Premium   -0.036764
L22.VIXCM30_Premium    0.041656
L23.VIXCM30_Premium   -0.023595
L24.VIXCM30_Premium   -0.014847
L25.VIXCM30_Premium    0.038173
L26.VIXCM30_Premium   -0.023949
L27.VIXCM30_Premium    0.057672
L28.VIXCM30_Premium   -0.028205
L29.VIXCM30_Premium   -0.055332
L30.VIXCM30_Premium    0.073213
L31.VIXCM30_Premium    0.022460
L32.VIXCM30_Premium   -0.031003
L33.VIXCM30_Premium   -0.053781
L34.VIXCM30_Premium   -0.021064
L35.VIXCM30_Premium    0.033376
L36.VIXCM30_Premium    0.013843
L37.VIXCM30_Premium    0.027811
L38.VIXCM30_Premium    0.023575
L39.VIXCM30_Premium   -0.042689
L40.VIXCM30_Premium   -0.013639
L41.VIXCM30_Premium   -0.025635
L42.VIXCM30_Premium    0.014546
L43.VIXCM30_Premium   -0.014341
L44.VIXCM30_Premium   -0.006181
L45.VIXCM30_Premium    0.014273
L46.VIXCM30_Premium   -0.003107
L47.VIXCM30_Premium    0.028452
L48.VIXCM30_Premium   -0.033419
L49.VIXCM30_Premium    0.034756
dtype: float64

print('RMSE (in sample):',np.sqrt(np.mean(model_fit.sigma2)))

RMSE (in sample): 1.1719350672913336

We will check to make sure that there is no information left in the residuals and that they are stationary:

run_ljung_box(model_fit.resid, lags=100)

Ljung-Box Test

plot_pacf(model_fit.resid, lags=40) # PACF looks good

model_fit.resid.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c1ce3e710>

run_ad_fuller(model_fit.resid)

ADF Statistic: -47.950249
p-value: 0.000000
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567

We reject the Null Hypothesis (H0) -- the time series is stationary

This model seems to caputure all the data so we will move forward.

Let's check to see how our model is performing on our test data:

H = train['VIXCM30_Premium'].tolist()
AR_ORDER = 49
AR_prem = []
AR_PRED=[]
VIX = []
ACTUAL_FUTURE=[]
ACTUAL_prem = []
n = 0
for d,row in train.tail(1).iterrows():
    VIX +=[row['VIX']]
for d,row in test.iterrows():
    ar_model = AR(H)
    ar_fit = ar_model.fit(maxlag=AR_ORDER)
    y_hat = ar_fit.predict(start=len(H), end=len(H) + 5)[0]
    AR_prem += [y_hat]
    AR_PRED+=[y_hat+VIX[-1]]
    VIX+=[row['VIX']]
    
    H += [row['VIXCM30_Premium']]
    ACTUAL_prem += [row['VIXCM30_Premium']]
    ACTUAL_FUTURE+=[row['VIXCM30']]
    n += 1

# plotting the first 200 points
plt.plot(AR_PRED[:200], label='predicted')
plt.plot(ACTUAL_FUTURE[:200], color='b', label='actual')
plt.legend()

<matplotlib.legend.Legend at 0x1c2e7a2da0>

plot_confusion_matrix(ACTUAL_FUTURE, AR_PRED)

Confusion Matrix
Raw
[[176 132]
 [160 120]]
Normalized
[[0.29931973 0.2244898 ]
 [0.27210884 0.20408163]]

{'cm': array([[176, 132],
        [160, 120]]), 'cm_norm': array([[0.29931973, 0.2244898 ],
        [0.27210884, 0.20408163]])}

ar_error2 = measure_error(ACTUAL_FUTURE, AR_PRED, label='AR(49)')
print(ar_error2)

RMSE = 0.8488649610470695
RMSE_pcent = 0.048114657402919224
label = AR(49)
{'RMSE': 0.8488649610470695, 'RMSE_pcent': 0.048114657402919224, 'label': 'AR(49)'}

Check to see how much money could be made out of sample:¶

make_money(ACTUAL_FUTURE, AR_PRED)

Total profit or loss: 5.9271
Number of correct days: 298
total incorrect days: 290

VIXCM60 Premium
¶

AR Model for VIXCM60 Premium¶

First we will look at the PACF plot of the premimums:

plot_pacf(train['VIXCM60_Premium'],lags=50)
plt.show()

# train autoregression
ar_model = AR(train['VIXCM60_Premium'])

Now we will check what order model is best from AIC and BIC scores:

aic = ar_model.select_order(maxlag=30,ic='aic',trend='nc')
bic = ar_model.select_order(maxlag=30,ic='bic',trend='nc')

print('model order (AIC): %d' % aic)
print('model order (BIC): %d' % bic)

model order (AIC): 30
model order (BIC): 5

Above we can see that the AIC score suggests 30 lags is the best, but this is the max lag we allowed it to use, so we will increase the max to see if it tells us something different.

aic = ar_model.select_order(maxlag=50,ic='aic',trend='nc')
bic = ar_model.select_order(maxlag=50,ic='bic',trend='nc')

print('model order (AIC): %d' % aic)
print('model order (BIC): %d' % bic)

model order (AIC): 35
model order (BIC): 5

We will try order =5 since it is the simplest model:

model_fit = ar_model.fit(maxlag=5)
print('Lag: %s' % model_fit.k_ar)
print('Coefficients: %s' % model_fit.params)
intercept = model_fit.params[0]
coeff = list(model_fit.params)[1:] # drop the intercept
coeff.reverse()

Lag: 5
Coefficients: const                 0.053017
L1.VIXCM60_Premium    0.769673
L2.VIXCM60_Premium    0.070797
L3.VIXCM60_Premium    0.064153
L4.VIXCM60_Premium   -0.024588
L5.VIXCM60_Premium    0.089369
dtype: float64

print('RMSE (in sample):',np.sqrt(np.mean(model_fit.sigma2)))

RMSE (in sample): 1.3704552438612214

We will check to make sure that there is no information left in the residuals and that they are stationary:

run_ljung_box(model_fit.resid, lags=40)

Ljung-Box Test
lag:  10 | p_value: 0.0001 | still has information
lag:  11 | p_value: 0.0002 | still has information
lag:  12 | p_value: 0.0002 | still has information
lag:  13 | p_value: 0.0004 | still has information
lag:  14 | p_value: 0.0002 | still has information
lag:  15 | p_value: 0.0004 | still has information
lag:  16 | p_value: 0.0002 | still has information
lag:  17 | p_value: 0.0003 | still has information
lag:  18 | p_value: 0.0001 | still has information
lag:  19 | p_value: 0.0001 | still has information
lag:  20 | p_value: 0.0002 | still has information
lag:  21 | p_value: 0.0000 | still has information
lag:  22 | p_value: 0.0001 | still has information
lag:  23 | p_value: 0.0001 | still has information
lag:  24 | p_value: 0.0001 | still has information
lag:  25 | p_value: 0.0002 | still has information
lag:  26 | p_value: 0.0002 | still has information
lag:  27 | p_value: 0.0001 | still has information
lag:  28 | p_value: 0.0001 | still has information
lag:  29 | p_value: 0.0000 | still has information
lag:  30 | p_value: 0.0000 | still has information
lag:  31 | p_value: 0.0000 | still has information
lag:  32 | p_value: 0.0000 | still has information
lag:  33 | p_value: 0.0000 | still has information
lag:  34 | p_value: 0.0000 | still has information
lag:  35 | p_value: 0.0000 | still has information
lag:  36 | p_value: 0.0000 | still has information
lag:  37 | p_value: 0.0000 | still has information
lag:  38 | p_value: 0.0000 | still has information
lag:  39 | p_value: 0.0000 | still has information
lag:  40 | p_value: 0.0000 | still has information

Now we will try order = 35 to see if we can caputure this missing information:

model_fit = ar_model.fit(maxlag=35)
print('Lag: %s' % model_fit.k_ar)
print('Coefficients: %s' % model_fit.params)
intercept = model_fit.params[0]
coeff = list(model_fit.params)[1:] # drop the intercept
coeff.reverse()

Lag: 35
Coefficients: const                  0.058063
L1.VIXCM60_Premium     0.770680
L2.VIXCM60_Premium     0.075459
L3.VIXCM60_Premium     0.058299
L4.VIXCM60_Premium    -0.020664
L5.VIXCM60_Premium     0.078253
L6.VIXCM60_Premium    -0.016905
L7.VIXCM60_Premium    -0.040313
L8.VIXCM60_Premium     0.073473
L9.VIXCM60_Premium     0.001452
L10.VIXCM60_Premium    0.092616
L11.VIXCM60_Premium   -0.074851
L12.VIXCM60_Premium   -0.030249
L13.VIXCM60_Premium   -0.005710
L14.VIXCM60_Premium    0.041038
L15.VIXCM60_Premium   -0.015750
L16.VIXCM60_Premium    0.033834
L17.VIXCM60_Premium   -0.022752
L18.VIXCM60_Premium   -0.054545
L19.VIXCM60_Premium    0.051477
L20.VIXCM60_Premium   -0.023131
L21.VIXCM60_Premium   -0.048491
L22.VIXCM60_Premium    0.049677
L23.VIXCM60_Premium   -0.005402
L24.VIXCM60_Premium   -0.010091
L25.VIXCM60_Premium    0.022448
L26.VIXCM60_Premium   -0.028361
L27.VIXCM60_Premium    0.051241
L28.VIXCM60_Premium   -0.018284
L29.VIXCM60_Premium   -0.052608
L30.VIXCM60_Premium    0.063715
L31.VIXCM60_Premium    0.017999
L32.VIXCM60_Premium   -0.035968
L33.VIXCM60_Premium   -0.036698
L34.VIXCM60_Premium   -0.025873
L35.VIXCM60_Premium    0.051801
dtype: float64

print('RMSE (in sample):',np.sqrt(np.mean(model_fit.sigma2)))

RMSE (in sample): 1.350898393202161

We will check to make sure that there is no information left in the residuals and that they are stationary:

run_ljung_box(model_fit.resid, lags=100)

Ljung-Box Test

plot_pacf(model_fit.resid, lags=40) # PACF looks good

model_fit.resid.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c1def65c0>

run_ad_fuller(model_fit.resid)

ADF Statistic: -48.165053
p-value: 0.000000
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567

We reject the Null Hypothesis (H0) -- the time series is stationary

This model seems to caputure all the data so we will move forward.

Let's check to see how our model is performing on our test data:

H = train['VIXCM60_Premium'].tolist()
AR_ORDER = 35
AR_prem1 = []
AR_PRED1=[]
VIX = []
ACTUAL_FUTURE1=[]
n = 0
for d,row in train.tail(1).iterrows():
    VIX +=[row['VIX']]
ACTUAL_prem1 = []
for d,row in test.iterrows():
    ar_model = AR(H)
    ar_fit = ar_model.fit(maxlag=AR_ORDER)
    y_hat = ar_fit.predict(start=len(H), end=len(H) + 5)[0]
    AR_prem1 += [y_hat]
    AR_PRED1+=[y_hat+VIX[-1]]
    VIX+=[row['VIX']]
    
    H += [row['VIXCM60_Premium']]
    ACTUAL_prem1 += [row['VIXCM60_Premium']]
    ACTUAL_FUTURE1+=[row['VIXCM60']]
    n += 1

# plotting the first 200 points
plt.plot(AR_PRED1[:200], label='predicted')
plt.plot(ACTUAL_FUTURE1[:200], color='b', label='actual')
plt.legend()

<matplotlib.legend.Legend at 0x1c1eec9a90>

plot_confusion_matrix(ACTUAL_FUTURE1, AR_PRED1)

Confusion Matrix
Raw
[[169 132]
 [165 122]]
Normalized
[[0.28741497 0.2244898 ]
 [0.28061224 0.20748299]]

{'cm': array([[169, 132],
        [165, 122]]), 'cm_norm': array([[0.28741497, 0.2244898 ],
        [0.28061224, 0.20748299]])}

ar_error2 = measure_error(ACTUAL_FUTURE1, AR_PRED1, label='AR(35)')
print(ar_error2)

RMSE = 0.6675266219320293
RMSE_pcent = 0.036184144108954146
label = AR(35)
{'RMSE': 0.6675266219320293, 'RMSE_pcent': 0.036184144108954146, 'label': 'AR(35)'}

Check to see how much money could be made out of sample:¶

make_money(ACTUAL_FUTURE1,AR_PRED1)

Total profit or loss: 0.9495
Number of correct days: 291
total incorrect days: 297

VIXCM90 Premium
¶

AR Model for VIXCM90 Premium¶

Now we will look at the PACF plot of the premimums:

plot_pacf(train['VIXCM90_Premium'],lags=50)
plt.show()

# train autoregression
ar_model = AR(train['VIXCM90_Premium'])

Now we will check what order model is best from AIC and BIC scores:

aic = ar_model.select_order(maxlag=30,ic='aic',trend='nc')
bic = ar_model.select_order(maxlag=30,ic='bic',trend='nc')

print('model order (AIC): %d' % aic)
print('model order (BIC): %d' % bic)

model order (AIC): 30
model order (BIC): 5

Above we can see that the AIC score suggests 30 lags is the best, but this is the max lag we allowed it to use, so we will increase the max to see if it tells us something different.

aic = ar_model.select_order(maxlag=50,ic='aic',trend='nc')
bic = ar_model.select_order(maxlag=50,ic='bic',trend='nc')

print('model order (AIC): %d' % aic)
print('model order (BIC): %d' % bic)

model order (AIC): 35
model order (BIC): 5

We will try order = 5 since it is the simplest model:

model_fit = ar_model.fit(maxlag=5)
print('Lag: %s' % model_fit.k_ar)
print('Coefficients: %s' % model_fit.params)
intercept = model_fit.params[0]
coeff = list(model_fit.params)[1:] # drop the intercept
coeff.reverse()

Lag: 5
Coefficients: const                 0.055619
L1.VIXCM90_Premium    0.769776
L2.VIXCM90_Premium    0.086111
L3.VIXCM90_Premium    0.058301
L4.VIXCM90_Premium   -0.030630
L5.VIXCM90_Premium    0.090508
dtype: float64

print('RMSE (in sample):',np.sqrt(np.mean(model_fit.sigma2)))

RMSE (in sample): 1.4572788232492344

We will check to make sure that there is no information left in the residuals and that they are stationary:

run_ljung_box(model_fit.resid, lags=40)

Ljung-Box Test
lag:  10 | p_value: 0.0001 | still has information
lag:  11 | p_value: 0.0001 | still has information
lag:  12 | p_value: 0.0001 | still has information
lag:  13 | p_value: 0.0002 | still has information
lag:  14 | p_value: 0.0002 | still has information
lag:  15 | p_value: 0.0003 | still has information
lag:  16 | p_value: 0.0001 | still has information
lag:  17 | p_value: 0.0001 | still has information
lag:  18 | p_value: 0.0001 | still has information
lag:  19 | p_value: 0.0001 | still has information
lag:  20 | p_value: 0.0001 | still has information
lag:  21 | p_value: 0.0000 | still has information
lag:  22 | p_value: 0.0001 | still has information
lag:  23 | p_value: 0.0001 | still has information
lag:  24 | p_value: 0.0001 | still has information
lag:  25 | p_value: 0.0001 | still has information
lag:  26 | p_value: 0.0002 | still has information
lag:  27 | p_value: 0.0001 | still has information
lag:  28 | p_value: 0.0001 | still has information
lag:  29 | p_value: 0.0000 | still has information
lag:  30 | p_value: 0.0000 | still has information
lag:  31 | p_value: 0.0000 | still has information
lag:  32 | p_value: 0.0000 | still has information
lag:  33 | p_value: 0.0000 | still has information
lag:  34 | p_value: 0.0000 | still has information
lag:  35 | p_value: 0.0000 | still has information
lag:  36 | p_value: 0.0000 | still has information
lag:  37 | p_value: 0.0000 | still has information
lag:  38 | p_value: 0.0000 | still has information
lag:  39 | p_value: 0.0000 | still has information
lag:  40 | p_value: 0.0000 | still has information

Now we will try order = 35 to see if we can caputure this missing information:

model_fit = ar_model.fit(maxlag=35)
print('Lag: %s' % model_fit.k_ar)
print('Coefficients: %s' % model_fit.params)
intercept = model_fit.params[0]
coeff = list(model_fit.params)[1:] # drop the intercept
coeff.reverse()

Lag: 35
Coefficients: const                  0.058823
L1.VIXCM90_Premium     0.768462
L2.VIXCM90_Premium     0.092559
L3.VIXCM90_Premium     0.053360
L4.VIXCM90_Premium    -0.026933
L5.VIXCM90_Premium     0.077405
L6.VIXCM90_Premium    -0.020206
L7.VIXCM90_Premium    -0.037089
L8.VIXCM90_Premium     0.069401
L9.VIXCM90_Premium     0.012646
L10.VIXCM90_Premium    0.089986
L11.VIXCM90_Premium   -0.076410
L12.VIXCM90_Premium   -0.035494
L13.VIXCM90_Premium   -0.006018
L14.VIXCM90_Premium    0.041574
L15.VIXCM90_Premium   -0.017509
L16.VIXCM90_Premium    0.042083
L17.VIXCM90_Premium   -0.021666
L18.VIXCM90_Premium   -0.061964
L19.VIXCM90_Premium    0.046902
L20.VIXCM90_Premium   -0.013907
L21.VIXCM90_Premium   -0.046341
L22.VIXCM90_Premium    0.040027
L23.VIXCM90_Premium    0.002175
L24.VIXCM90_Premium   -0.013742
L25.VIXCM90_Premium    0.025134
L26.VIXCM90_Premium   -0.027931
L27.VIXCM90_Premium    0.050044
L28.VIXCM90_Premium   -0.010551
L29.VIXCM90_Premium   -0.052993
L30.VIXCM90_Premium    0.061218
L31.VIXCM90_Premium    0.016391
L32.VIXCM90_Premium   -0.040458
L33.VIXCM90_Premium   -0.034215
L34.VIXCM90_Premium   -0.024925
L35.VIXCM90_Premium    0.051781
dtype: float64

print('RMSE (in sample):',np.sqrt(np.mean(model_fit.sigma2)))

RMSE (in sample): 1.4367561404498372

We will check to make sure that there is no information left in the residuals and that they are stationary:

run_ljung_box(model_fit.resid, lags=100)

Ljung-Box Test

plot_pacf(model_fit.resid, lags=40) # PACF looks good

model_fit.resid.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c2c267f98>

run_ad_fuller(model_fit.resid)

ADF Statistic: -48.148136
p-value: 0.000000
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567

We reject the Null Hypothesis (H0) -- the time series is stationary

This model seems to caputure all the data so we will move forward.

Let's check to see how our model is performing on our test data:

H = train['VIXCM90_Premium'].tolist()
AR_ORDER = 35
AR_prem2 = []
AR_PRED2=[]
VIX = []
ACTUAL_FUTURE2=[]
n = 0
for d,row in train.tail(1).iterrows():
    VIX +=[row['VIX']]
ACTUAL_prem2 = []
for d,row in test.iterrows():
    ar_model = AR(H)
    ar_fit = ar_model.fit(maxlag=AR_ORDER)
    y_hat = ar_fit.predict(start=len(H), end=len(H) + 5)[0]
    AR_prem2 += [y_hat]
    AR_PRED2+=[y_hat+VIX[-1]]
    VIX+=[row['VIX']]
    
    H += [row['VIXCM90_Premium']]
    ACTUAL_prem2 += [row['VIXCM90_Premium']]
    ACTUAL_FUTURE2+=[row['VIXCM90']]
    n += 1

# plotting the first 200 points
plt.plot(AR_PRED2[:200], label='predicted')
plt.plot(ACTUAL_FUTURE2[:200], color='b', label='actual')
plt.legend()

<matplotlib.legend.Legend at 0x1c2c302fd0>

plot_confusion_matrix(ACTUAL_FUTURE2, AR_PRED2)

Confusion Matrix
Raw
[[160 134]
 [168 126]]
Normalized
[[0.27210884 0.22789116]
 [0.28571429 0.21428571]]

{'cm': array([[160, 134],
        [168, 126]]), 'cm_norm': array([[0.27210884, 0.22789116],
        [0.28571429, 0.21428571]])}

ar_error2 = measure_error(ACTUAL_FUTURE2, AR_PRED2, label='AR(35)')
print(ar_error2)

RMSE = 0.6043100435289301
RMSE_pcent = 0.031809285219358684
label = AR(35)
{'RMSE': 0.6043100435289301, 'RMSE_pcent': 0.031809285219358684, 'label': 'AR(35)'}

Check to see how much money could be made out of sample:¶

make_money(ACTUAL_FUTURE2, AR_PRED2)

Total profit or loss: 2.9856
Number of correct days: 286
total incorrect days: 302

VIXCM120 Premium
¶

VIXCM120 Premium¶

First we will look at the PACF plot of the premimums:

plot_pacf(train['VIXCM120_Premium'],lags=50)
plt.show()

# train autoregression
ar_model = AR(train['VIXCM120_Premium'])

Now we will check what order model is best from AIC and BIC scores:

aic = ar_model.select_order(maxlag=30,ic='aic',trend='nc')
bic = ar_model.select_order(maxlag=30,ic='bic',trend='nc')

print('model order (AIC): %d' % aic)
print('model order (BIC): %d' % bic)

model order (AIC): 30
model order (BIC): 5

Above we can see that the AIC score suggests 30 lags is the best, but this is the max lag we allowed it to use, so we will increase the max to see if it tells us something different.

aic = ar_model.select_order(maxlag=50,ic='aic',trend='nc')
bic = ar_model.select_order(maxlag=50,ic='bic',trend='nc')

print('model order (AIC): %d' % aic)
print('model order (BIC): %d' % bic)

model order (AIC): 35
model order (BIC): 5

We will try order = 5 since it is the simplest model:

model_fit = ar_model.fit(maxlag=5)
print('Lag: %s' % model_fit.k_ar)
print('Coefficients: %s' % model_fit.params)
intercept = model_fit.params[0]
coeff = list(model_fit.params)[1:] # drop the intercept
coeff.reverse()

Lag: 5
Coefficients: const                  0.058267
L1.VIXCM120_Premium    0.766850
L2.VIXCM120_Premium    0.089744
L3.VIXCM120_Premium    0.064586
L4.VIXCM120_Premium   -0.040856
L5.VIXCM120_Premium    0.096360
dtype: float64

print('RMSE (in sample):',np.sqrt(np.mean(model_fit.sigma2)))

RMSE (in sample): 1.5162838219095407

We will check to make sure that there is no information left in the residuals and that they are stationary:

run_ljung_box(model_fit.resid, lags=100)

Ljung-Box Test
lag:  10 | p_value: 0.0000 | still has information
lag:  11 | p_value: 0.0000 | still has information
lag:  12 | p_value: 0.0000 | still has information
lag:  13 | p_value: 0.0001 | still has information
lag:  14 | p_value: 0.0001 | still has information
lag:  15 | p_value: 0.0001 | still has information
lag:  16 | p_value: 0.0000 | still has information
lag:  17 | p_value: 0.0000 | still has information
lag:  18 | p_value: 0.0000 | still has information
lag:  19 | p_value: 0.0000 | still has information
lag:  20 | p_value: 0.0000 | still has information
lag:  21 | p_value: 0.0000 | still has information
lag:  22 | p_value: 0.0000 | still has information
lag:  23 | p_value: 0.0000 | still has information
lag:  24 | p_value: 0.0000 | still has information
lag:  25 | p_value: 0.0000 | still has information
lag:  26 | p_value: 0.0000 | still has information
lag:  27 | p_value: 0.0000 | still has information
lag:  28 | p_value: 0.0000 | still has information
lag:  29 | p_value: 0.0000 | still has information
lag:  30 | p_value: 0.0000 | still has information
lag:  31 | p_value: 0.0000 | still has information
lag:  32 | p_value: 0.0000 | still has information
lag:  33 | p_value: 0.0000 | still has information
lag:  34 | p_value: 0.0000 | still has information
lag:  35 | p_value: 0.0000 | still has information
lag:  36 | p_value: 0.0000 | still has information
lag:  37 | p_value: 0.0000 | still has information
lag:  38 | p_value: 0.0000 | still has information
lag:  39 | p_value: 0.0000 | still has information
lag:  40 | p_value: 0.0000 | still has information
lag:  41 | p_value: 0.0000 | still has information
lag:  42 | p_value: 0.0000 | still has information
lag:  43 | p_value: 0.0000 | still has information
lag:  44 | p_value: 0.0000 | still has information
lag:  45 | p_value: 0.0000 | still has information
lag:  46 | p_value: 0.0000 | still has information
lag:  47 | p_value: 0.0000 | still has information
lag:  48 | p_value: 0.0000 | still has information
lag:  49 | p_value: 0.0000 | still has information
lag:  50 | p_value: 0.0000 | still has information
lag:  51 | p_value: 0.0000 | still has information
lag:  52 | p_value: 0.0000 | still has information
lag:  53 | p_value: 0.0000 | still has information
lag:  54 | p_value: 0.0000 | still has information
lag:  55 | p_value: 0.0000 | still has information
lag:  56 | p_value: 0.0000 | still has information
lag:  57 | p_value: 0.0000 | still has information
lag:  58 | p_value: 0.0000 | still has information
lag:  59 | p_value: 0.0000 | still has information
lag:  60 | p_value: 0.0000 | still has information
lag:  61 | p_value: 0.0000 | still has information
lag:  62 | p_value: 0.0000 | still has information
lag:  63 | p_value: 0.0000 | still has information
lag:  64 | p_value: 0.0000 | still has information
lag:  65 | p_value: 0.0000 | still has information
lag:  66 | p_value: 0.0000 | still has information
lag:  67 | p_value: 0.0000 | still has information
lag:  68 | p_value: 0.0000 | still has information
lag:  69 | p_value: 0.0000 | still has information
lag:  70 | p_value: 0.0000 | still has information
lag:  71 | p_value: 0.0000 | still has information
lag:  72 | p_value: 0.0000 | still has information
lag:  73 | p_value: 0.0000 | still has information
lag:  74 | p_value: 0.0000 | still has information
lag:  75 | p_value: 0.0000 | still has information
lag:  76 | p_value: 0.0000 | still has information
lag:  77 | p_value: 0.0000 | still has information
lag:  78 | p_value: 0.0000 | still has information
lag:  79 | p_value: 0.0000 | still has information
lag:  80 | p_value: 0.0000 | still has information
lag:  81 | p_value: 0.0000 | still has information
lag:  82 | p_value: 0.0000 | still has information
lag:  83 | p_value: 0.0000 | still has information
lag:  84 | p_value: 0.0000 | still has information
lag:  85 | p_value: 0.0000 | still has information
lag:  86 | p_value: 0.0000 | still has information
lag:  87 | p_value: 0.0000 | still has information
lag:  88 | p_value: 0.0000 | still has information
lag:  89 | p_value: 0.0000 | still has information
lag:  90 | p_value: 0.0000 | still has information
lag:  91 | p_value: 0.0000 | still has information
lag:  92 | p_value: 0.0000 | still has information
lag:  93 | p_value: 0.0000 | still has information
lag:  94 | p_value: 0.0000 | still has information
lag:  95 | p_value: 0.0000 | still has information
lag:  96 | p_value: 0.0000 | still has information
lag:  97 | p_value: 0.0000 | still has information
lag:  98 | p_value: 0.0000 | still has information
lag:  99 | p_value: 0.0000 | still has information
lag: 100 | p_value: 0.0000 | still has information

Now we will try order = 35 to see if we can caputure this missing information:

model_fit = ar_model.fit(maxlag=35)
print('Lag: %s' % model_fit.k_ar)
print('Coefficients: %s' % model_fit.params)
intercept = model_fit.params[0]
coeff = list(model_fit.params)[1:] # drop the intercept
coeff.reverse()

Lag: 35
Coefficients: const                   0.062413
L1.VIXCM120_Premium     0.768963
L2.VIXCM120_Premium     0.095499
L3.VIXCM120_Premium     0.056335
L4.VIXCM120_Premium    -0.033309
L5.VIXCM120_Premium     0.083427
L6.VIXCM120_Premium    -0.015385
L7.VIXCM120_Premium    -0.042856
L8.VIXCM120_Premium     0.070788
L9.VIXCM120_Premium     0.005830
L10.VIXCM120_Premium    0.099565
L11.VIXCM120_Premium   -0.087578
L12.VIXCM120_Premium   -0.040525
L13.VIXCM120_Premium    0.003309
L14.VIXCM120_Premium    0.032214
L15.VIXCM120_Premium   -0.019641
L16.VIXCM120_Premium    0.051833
L17.VIXCM120_Premium   -0.025325
L18.VIXCM120_Premium   -0.063755
L19.VIXCM120_Premium    0.058128
L20.VIXCM120_Premium   -0.012063
L21.VIXCM120_Premium   -0.054492
L22.VIXCM120_Premium    0.052034
L23.VIXCM120_Premium   -0.002868
L24.VIXCM120_Premium   -0.017726
L25.VIXCM120_Premium    0.027768
L26.VIXCM120_Premium   -0.022699
L27.VIXCM120_Premium    0.039281
L28.VIXCM120_Premium   -0.000798
L29.VIXCM120_Premium   -0.058814
L30.VIXCM120_Premium    0.061087
L31.VIXCM120_Premium    0.009791
L32.VIXCM120_Premium   -0.033172
L33.VIXCM120_Premium   -0.034713
L34.VIXCM120_Premium   -0.024950
L35.VIXCM120_Premium    0.049975
dtype: float64

print('RMSE (in sample):',np.sqrt(np.mean(model_fit.sigma2)))

RMSE (in sample): 1.4936875650982426

We will check to make sure that there is no information left in the residuals and that they are stationary:

run_ljung_box(model_fit.resid, lags=100)

Ljung-Box Test

plot_pacf(model_fit.resid, lags=40) # PACF looks good

model_fit.resid.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c2e3a7470>

run_ad_fuller(model_fit.resid)

ADF Statistic: -48.151979
p-value: 0.000000
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567

We reject the Null Hypothesis (H0) -- the time series is stationary

This model seems to caputure all the data so we will move forward.

Let's check to see how our model is performing on our test data:

H = train['VIXCM120_Premium'].tolist()
AR_ORDER = 35
AR_prem3 = []
AR_PRED3=[]
VIX = []
ACTUAL_FUTURE3=[]
n = 0
for d,row in train.tail(1).iterrows():
    VIX +=[row['VIX']]
ACTUAL_prem3 = []
for d,row in test.iterrows():
    ar_model = AR(H)
    ar_fit = ar_model.fit(maxlag=AR_ORDER)
    y_hat = ar_fit.predict(start=len(H), end=len(H) + 5)[0]
    AR_prem3 += [y_hat]
    AR_PRED3+=[y_hat+VIX[-1]]
    VIX+=[row['VIX']]
    
    H += [row['VIXCM120_Premium']]
    ACTUAL_prem3 += [row['VIXCM120_Premium']]
    ACTUAL_FUTURE3+=[row['VIXCM120']]
    n += 1

# plotting the first 200 points
plt.plot(AR_PRED3[:200], label='predicted')
plt.plot(ACTUAL_FUTURE3[:200], color='b', label='actual')
plt.legend()

<matplotlib.legend.Legend at 0x1c1d24e588>

plot_confusion_matrix(ACTUAL_FUTURE3, AR_PRED3)

Confusion Matrix
Raw
[[154 132]
 [173 129]]
Normalized
[[0.26190476 0.2244898 ]
 [0.29421769 0.21938776]]

{'cm': array([[154, 132],
        [173, 129]]), 'cm_norm': array([[0.26190476, 0.2244898 ],
        [0.29421769, 0.21938776]])}

ar_error2 = measure_error(ACTUAL_FUTURE3, AR_PRED3, label='AR(35)')
print(ar_error2)

RMSE = 0.5713382942295682
RMSE_pcent = 0.029363389373415195
label = AR(35)
{'RMSE': 0.5713382942295682, 'RMSE_pcent': 0.029363389373415195, 'label': 'AR(35)'}

Check to see how much money could be made out of sample:¶

make_money(ACTUAL_FUTURE3, AR_PRED3)

Total profit or loss: 4.7723
Number of correct days: 283
total incorrect days: 305

VIXCM150 Premium
¶

VIXCM150 Premium¶

First we will look at the PACF plot of the premimums:

plot_pacf(train['VIXCM150_Premium'],lags=50)
plt.show()

# train autoregression
ar_model = AR(train['VIXCM150_Premium'])

Now we will check what order model is best from AIC and BIC scores:

aic = ar_model.select_order(maxlag=30,ic='aic',trend='nc')
bic = ar_model.select_order(maxlag=30,ic='bic',trend='nc')

print('model order (AIC): %d' % aic)
print('model order (BIC): %d' % bic)

model order (AIC): 30
model order (BIC): 5

Above we can see that the AIC score suggests 30 lags is the best, but this is the max lag we allowed it to use, so we will increase the max to see if it tells us something different.

aic = ar_model.select_order(maxlag=50,ic='aic',trend='nc')
bic = ar_model.select_order(maxlag=50,ic='bic',trend='nc')

print('model order (AIC): %d' % aic)
print('model order (BIC): %d' % bic)

model order (AIC): 49
model order (BIC): 11

We will try order = 11 since it is the simplest model:

model_fit = ar_model.fit(maxlag=11)
print('Lag: %s' % model_fit.k_ar)
print('Coefficients: %s' % model_fit.params)
intercept = model_fit.params[0]
coeff = list(model_fit.params)[1:] # drop the intercept
coeff.reverse()

Lag: 11
Coefficients: const                   0.061096
L1.VIXCM150_Premium     0.767740
L2.VIXCM150_Premium     0.091239
L3.VIXCM150_Premium     0.072728
L4.VIXCM150_Premium    -0.054116
L5.VIXCM150_Premium     0.098360
L6.VIXCM150_Premium    -0.000875
L7.VIXCM150_Premium    -0.067692
L8.VIXCM150_Premium     0.080361
L9.VIXCM150_Premium    -0.000225
L10.VIXCM150_Premium    0.099617
L11.VIXCM150_Premium   -0.108832
dtype: float64

print('RMSE (in sample):',np.sqrt(np.mean(model_fit.sigma2)))

RMSE (in sample): 1.5721360055497136

We will check to make sure that there is no information left in the residuals and that they are stationary:

run_ljung_box(model_fit.resid, lags=100)

Ljung-Box Test
lag:  30 | p_value: 0.0346 | still has information
lag:  31 | p_value: 0.0362 | still has information
lag:  32 | p_value: 0.0466 | still has information
lag:  34 | p_value: 0.0090 | still has information
lag:  35 | p_value: 0.0120 | still has information
lag:  36 | p_value: 0.0143 | still has information
lag:  37 | p_value: 0.0186 | still has information
lag:  38 | p_value: 0.0093 | still has information
lag:  39 | p_value: 0.0101 | still has information
lag:  40 | p_value: 0.0130 | still has information
lag:  41 | p_value: 0.0155 | still has information
lag:  42 | p_value: 0.0184 | still has information
lag:  43 | p_value: 0.0092 | still has information
lag:  44 | p_value: 0.0118 | still has information
lag:  45 | p_value: 0.0120 | still has information
lag:  46 | p_value: 0.0142 | still has information
lag:  47 | p_value: 0.0132 | still has information
lag:  48 | p_value: 0.0010 | still has information
lag:  49 | p_value: 0.0014 | still has information
lag:  50 | p_value: 0.0001 | still has information
lag:  51 | p_value: 0.0001 | still has information
lag:  52 | p_value: 0.0001 | still has information
lag:  53 | p_value: 0.0001 | still has information
lag:  54 | p_value: 0.0001 | still has information
lag:  55 | p_value: 0.0001 | still has information
lag:  56 | p_value: 0.0000 | still has information
lag:  57 | p_value: 0.0000 | still has information
lag:  58 | p_value: 0.0000 | still has information
lag:  59 | p_value: 0.0000 | still has information
lag:  60 | p_value: 0.0000 | still has information
lag:  61 | p_value: 0.0000 | still has information
lag:  62 | p_value: 0.0000 | still has information
lag:  63 | p_value: 0.0000 | still has information
lag:  64 | p_value: 0.0000 | still has information
lag:  65 | p_value: 0.0000 | still has information
lag:  66 | p_value: 0.0000 | still has information
lag:  67 | p_value: 0.0000 | still has information
lag:  68 | p_value: 0.0000 | still has information
lag:  69 | p_value: 0.0000 | still has information
lag:  70 | p_value: 0.0000 | still has information
lag:  71 | p_value: 0.0000 | still has information
lag:  72 | p_value: 0.0000 | still has information
lag:  73 | p_value: 0.0000 | still has information
lag:  74 | p_value: 0.0000 | still has information
lag:  75 | p_value: 0.0000 | still has information
lag:  76 | p_value: 0.0000 | still has information
lag:  77 | p_value: 0.0000 | still has information
lag:  78 | p_value: 0.0000 | still has information
lag:  79 | p_value: 0.0000 | still has information
lag:  80 | p_value: 0.0000 | still has information
lag:  81 | p_value: 0.0000 | still has information
lag:  82 | p_value: 0.0000 | still has information
lag:  83 | p_value: 0.0000 | still has information
lag:  84 | p_value: 0.0000 | still has information
lag:  85 | p_value: 0.0000 | still has information
lag:  86 | p_value: 0.0000 | still has information
lag:  87 | p_value: 0.0000 | still has information
lag:  88 | p_value: 0.0000 | still has information
lag:  89 | p_value: 0.0000 | still has information
lag:  90 | p_value: 0.0000 | still has information
lag:  91 | p_value: 0.0000 | still has information
lag:  92 | p_value: 0.0000 | still has information
lag:  93 | p_value: 0.0000 | still has information
lag:  94 | p_value: 0.0000 | still has information
lag:  95 | p_value: 0.0000 | still has information
lag:  96 | p_value: 0.0000 | still has information
lag:  97 | p_value: 0.0000 | still has information
lag:  98 | p_value: 0.0000 | still has information
lag:  99 | p_value: 0.0000 | still has information
lag: 100 | p_value: 0.0000 | still has information

Now we will try order = 49 to see if we can caputure this missing information:

model_fit = ar_model.fit(maxlag=49)
print('Lag: %s' % model_fit.k_ar)
print('Coefficients: %s' % model_fit.params)
intercept = model_fit.params[0]
coeff = list(model_fit.params)[1:] # drop the intercept
coeff.reverse()

Lag: 49
Coefficients: const                   0.061411
L1.VIXCM150_Premium     0.769016
L2.VIXCM150_Premium     0.096952
L3.VIXCM150_Premium     0.062871
L4.VIXCM150_Premium    -0.037810
L5.VIXCM150_Premium     0.087700
L6.VIXCM150_Premium    -0.010294
L7.VIXCM150_Premium    -0.053452
L8.VIXCM150_Premium     0.072443
L9.VIXCM150_Premium     0.000395
L10.VIXCM150_Premium    0.108985
L11.VIXCM150_Premium   -0.098503
L12.VIXCM150_Premium   -0.045901
L13.VIXCM150_Premium    0.011974
L14.VIXCM150_Premium    0.023194
L15.VIXCM150_Premium   -0.017758
L16.VIXCM150_Premium    0.063046
L17.VIXCM150_Premium   -0.032031
L18.VIXCM150_Premium   -0.062360
L19.VIXCM150_Premium    0.062948
L20.VIXCM150_Premium   -0.006435
L21.VIXCM150_Premium   -0.060082
L22.VIXCM150_Premium    0.056023
L23.VIXCM150_Premium   -0.002980
L24.VIXCM150_Premium   -0.024719
L25.VIXCM150_Premium    0.030095
L26.VIXCM150_Premium   -0.012988
L27.VIXCM150_Premium    0.026033
L28.VIXCM150_Premium    0.001659
L29.VIXCM150_Premium   -0.057286
L30.VIXCM150_Premium    0.056237
L31.VIXCM150_Premium    0.006627
L32.VIXCM150_Premium   -0.021767
L33.VIXCM150_Premium   -0.036470
L34.VIXCM150_Premium   -0.026932
L35.VIXCM150_Premium    0.039700
L36.VIXCM150_Premium    0.006839
L37.VIXCM150_Premium    0.016915
L38.VIXCM150_Premium    0.027289
L39.VIXCM150_Premium   -0.011537
L40.VIXCM150_Premium   -0.014068
L41.VIXCM150_Premium   -0.034881
L42.VIXCM150_Premium    0.004053
L43.VIXCM150_Premium   -0.013500
L44.VIXCM150_Premium    0.016294
L45.VIXCM150_Premium   -0.006829
L46.VIXCM150_Premium   -0.006763
L47.VIXCM150_Premium    0.034234
L48.VIXCM150_Premium   -0.063235
L49.VIXCM150_Premium    0.055298
dtype: float64

print('RMSE (in sample):',np.sqrt(np.mean(model_fit.sigma2)))

RMSE (in sample): 1.5569022259180805

We will check to make sure that there is no information left in the residuals and that they are stationary:

run_ljung_box(model_fit.resid, lags=100)

Ljung-Box Test

plot_pacf(model_fit.resid, lags=40) # PACF looks good

model_fit.resid.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c2e27ee48>

run_ad_fuller(model_fit.resid)

ADF Statistic: -47.991732
p-value: 0.000000
Critical Values:
	1%: -3.433
	5%: -2.863
	10%: -2.567

We reject the Null Hypothesis (H0) -- the time series is stationary

This model seems to caputure all the data so we will move forward.

Let's check to see how our model is performing on our test data:

H = train['VIXCM150_Premium'].tolist()
AR_ORDER = 49
AR_prem4 = []
AR_PRED4=[]
VIX = []
ACTUAL_FUTURE4=[]
n = 0
for d,row in train.tail(1).iterrows():
    VIX +=[row['VIX']]
ACTUAL_prem4 = []
for d,row in test.iterrows():
    ar_model = AR(H)
    ar_fit = ar_model.fit(maxlag=AR_ORDER)
    y_hat = ar_fit.predict(start=len(H), end=len(H) + 5)[0]
    AR_prem4 += [y_hat]
    AR_PRED4+=[y_hat+VIX[-1]]
    VIX+=[row['VIX']]
    
    H += [row['VIXCM150_Premium']]
    ACTUAL_prem4 += [row['VIXCM150_Premium']]
    ACTUAL_FUTURE4+=[row['VIXCM150']]
    n += 1

# plotting the first 200 points
plt.plot(AR_PRED4[:200], label='predicted')
plt.plot(ACTUAL_FUTURE4[:200], color='b', label='actual')
plt.legend()

<matplotlib.legend.Legend at 0x1c2ec1b7f0>

plot_confusion_matrix(ACTUAL_FUTURE4, AR_PRED4)

Confusion Matrix
Raw
[[153 135]
 [159 141]]
Normalized
[[0.26020408 0.22959184]
 [0.27040816 0.23979592]]

{'cm': array([[153, 135],
        [159, 141]]), 'cm_norm': array([[0.26020408, 0.22959184],
        [0.27040816, 0.23979592]])}

ar_error2 = measure_error(ACTUAL_FUTURE4, AR_PRED4, label='AR(49)')
print(ar_error2)

RMSE = 0.5579351162667225
RMSE_pcent = 0.02801621318220747
label = AR(49)
{'RMSE': 0.5579351162667225, 'RMSE_pcent': 0.02801621318220747, 'label': 'AR(49)'}

Check to see how much money could be made out of sample:¶

make_money(ACTUAL_FUTURE4, AR_PRED4)

Total profit or loss: 5.1826
Number of correct days: 294
total incorrect days: 294

Summary of AR Models
¶

Summary of AR Models:¶

data	model	RMSE_Percent	PnL
VIXCM30_Premium	AR(49)	0.04811	+\$ 5.9271
VIXCM60_Premium	AR(35)	0.03618	+\$ 0.9495
VIXCM90_Premium	AR(35)	0.03180	+\$ 2.9856
VIXCM120_Premium	AR(35)	0.02930	+\$ 4.7723
VIXCM150_Premium	AR(49)	0.02801	+\$ 5.1826

Note - we did try to add rolling forecasts, limiting the look back to 180, 90, and 60 day, but this actually made the model perform worse, so we decided to keep all historical data.

ARMA
¶

ARMA Models¶

*Note we will conclude that all variables used in the following model are stationary as we already checked above.

ARMA Model for VIXCM30 Premium¶

First let's take a look at the PACF and ACF:

X= train['VIXCM30_Premium']

auto = acf(X)
print(auto)
y = plot_acf(X, lags=60)
plt.show()
y = plot_pacf(X, lags=40)
plt.show()

[1.         0.93263131 0.89195702 0.86420381 0.84062773 0.82582195
 0.80192867 0.7817517  0.77243611 0.75746196 0.74231963 0.71630187
 0.68909109 0.66980097 0.65316705 0.63086905 0.60978293 0.58550488
 0.55986135 0.54301903 0.519158   0.49584697 0.47905112 0.45706396
 0.43962203 0.42561678 0.40758918 0.39513983 0.37418473 0.35365751
 0.34368655 0.32654458 0.30575451 0.28431679 0.26717804 0.26069954
 0.25273016 0.24239629 0.23019658 0.21311427 0.1994161 ]

We will runa grid search to help us decide what p and q values to use in our model:

grid_result = run_grid_over_arma(X, p_lags=list(range(5)), q_lags=range(5))

['p(0)', 'n/a', '10572.12', ' 9608.72', ' 8927.59', ' 8681.60']
['p(1)', ' 7682.91', ' 7595.44', ' 7570.22', ' 7573.93', ' 7576.89']
['p(2)', ' 7619.48', ' 7564.96', ' 7572.30', ' 7580.01', 'n/a']
['p(3)', ' 7591.93', ' 7572.28', ' 7578.77', ' 7586.41', ' 7565.84']
['p(4)', ' 7588.03', ' 7579.89', ' 7587.41', 'n/a', 'n/a']
+------+----------+----------+----------+----------+----------+
| p\q  |  q (0)   |  q (1)   |  q (2)   |  q (3)   |  q (4)   |
+------+----------+----------+----------+----------+----------+
| p(0) |   n/a    | 10572.12 |  9608.72 |  8927.59 |  8681.60 |
| p(1) |  7682.91 |  7595.44 |  7570.22 |  7573.93 |  7576.89 |
| p(2) |  7619.48 |  7564.96 |  7572.30 |  7580.01 |   n/a    |
| p(3) |  7591.93 |  7572.28 |  7578.77 |  7586.41 |  7565.84 |
| p(4) |  7588.03 |  7579.89 |  7587.41 |   n/a    |   n/a    |
+------+----------+----------+----------+----------+----------+
bad values: [(0, 0), (2, 4), (4, 3), (4, 4)]
time required 9.743584871292114

b = sorted(list(grid_result.items()), key=lambda k: k[1])
print(b)

[((2, 1), 7564.957472965394), ((3, 4), 7565.8438456627155), ((1, 2), 7570.21967803622), ((3, 1), 7572.279762769047), ((2, 2), 7572.299550456719), ((1, 3), 7573.930732382952), ((1, 4), 7576.894075550767), ((3, 2), 7578.7665201778), ((4, 1), 7579.892549913009), ((2, 3), 7580.011855675479), ((3, 3), 7586.414384936876), ((4, 2), 7587.408023579574), ((4, 0), 7588.029367390866), ((3, 0), 7591.9328531265255), ((1, 1), 7595.437259910496), ((2, 0), 7619.4787994317485), ((1, 0), 7682.908739805367), ((0, 4), 8681.60204393116), ((0, 3), 8927.590315352254), ((0, 2), 9608.724979676586), ((0, 1), 10572.123910239005)]

We will try an ARMA(2,1) model and a ARMA (3,4) model:

First ARMA(2,1):

model = ARIMA(X, order=(2,0,1))
r = model.fit()
r.summary()
print('RMSE (in sample):',np.sqrt(np.mean(r.sigma2)))

RMSE (in sample): 1.1975193171875995

train1=train['VIXCM30_Premium'].dropna().tolist()
arima_model = ARIMA(train1, order=(2,0,1))
model_fit = arima_model.fit()

We can check to see if there is any information left in the residuals:

run_ljung_box(model_fit.resid,lags=40)

Ljung-Box Test
lag:   5 | p_value: 0.0245 | still has information
lag:   6 | p_value: 0.0342 | still has information
lag:   7 | p_value: 0.0024 | still has information
lag:   8 | p_value: 0.0006 | still has information
lag:   9 | p_value: 0.0005 | still has information
lag:  10 | p_value: 0.0000 | still has information
lag:  11 | p_value: 0.0000 | still has information
lag:  12 | p_value: 0.0000 | still has information
lag:  13 | p_value: 0.0000 | still has information
lag:  14 | p_value: 0.0000 | still has information
lag:  15 | p_value: 0.0000 | still has information
lag:  16 | p_value: 0.0000 | still has information
lag:  17 | p_value: 0.0000 | still has information
lag:  18 | p_value: 0.0000 | still has information
lag:  19 | p_value: 0.0000 | still has information
lag:  20 | p_value: 0.0000 | still has information
lag:  21 | p_value: 0.0000 | still has information
lag:  22 | p_value: 0.0000 | still has information
lag:  23 | p_value: 0.0000 | still has information
lag:  24 | p_value: 0.0000 | still has information
lag:  25 | p_value: 0.0000 | still has information
lag:  26 | p_value: 0.0000 | still has information
lag:  27 | p_value: 0.0000 | still has information
lag:  28 | p_value: 0.0000 | still has information
lag:  29 | p_value: 0.0000 | still has information
lag:  30 | p_value: 0.0000 | still has information
lag:  31 | p_value: 0.0000 | still has information
lag:  32 | p_value: 0.0000 | still has information
lag:  33 | p_value: 0.0000 | still has information
lag:  34 | p_value: 0.0000 | still has information
lag:  35 | p_value: 0.0000 | still has information
lag:  36 | p_value: 0.0000 | still has information
lag:  37 | p_value: 0.0000 | still has information
lag:  38 | p_value: 0.0000 | still has information
lag:  39 | p_value: 0.0000 | still has information
lag:  40 | p_value: 0.0000 | still has information

Unfortunately, there is still a lot of data left in the residuals. We will check to see how this model performs out of sample regardless and then we will check to see if a model with a different order might capture more of this residual information.

H = train['VIXCM30_Premium'].tolist()
ARIMA_ORDER = (2,0,1)
ARMA_prem4 = []
ARMA_PRED4=[]
VIX = []
ACTUAL_FUTURE4=[]
ACTUAL_prem4 = []
n = 0

for d,row in train.tail(1).iterrows():
    VIX +=[row['VIX']]

for d,row in test.iterrows():
    arima_model = ARIMA(H, order=ARIMA_ORDER)
    arima_fit = arima_model.fit()
    y_hat = arima_fit.predict(start=len(H), end=len(H) + 5)[0]
    ARMA_prem4 += [y_hat]
    ARMA_PRED4+=[y_hat+VIX[-1]]
    VIX+=[row['VIX']]
    
    H += [row['VIXCM30_Premium']]
    ACTUAL_prem4 += [row['VIXCM30_Premium']]
    ACTUAL_FUTURE4+=[row['VIXCM30']]
    n += 1


# plotting the first 200 points
plt.plot(ARMA_PRED4[:200], label='predicted')
plt.plot(ACTUAL_FUTURE4[:200], color='b', label='actual')
plt.legend()

<matplotlib.legend.Legend at 0x1c2f5335c0>

rolling_ar_errors = measure_error(ACTUAL_FUTURE4, ARMA_PRED4, label='Rolling ARIMA(%d,%d,%d)' % ARIMA_ORDER)

RMSE = 0.8369772457435734
RMSE_pcent = 0.047440847815554896
label = Rolling ARIMA(2,0,1)

ar_cm = plot_confusion_matrix(ACTUAL_FUTURE4, ARMA_PRED4)

Confusion Matrix
Raw
[[193 115]
 [181  99]]
Normalized
[[0.32823129 0.19557823]
 [0.30782313 0.16836735]]

make_money(ACTUAL_FUTURE4, ARMA_PRED4)

Total profit or loss: 1.5630
Number of correct days: 294
total incorrect days: 294

Now let's try our ARMA (3,4) model to see if we can solve the issue of the information that is remaining in the residuals:

model2 = ARIMA(X, order=(3,0,4))
r2 = model2.fit()
r2.summary()
print('RMSE (in sample):',np.sqrt(np.mean(r2.sigma2)))

RMSE (in sample): 1.1897402525654754

arima_model = ARIMA(train1, order=(3,0,4))
model_fit2 = arima_model.fit()

run_ljung_box(r2.resid,lags=40)

Ljung-Box Test
lag:  10 | p_value: 0.0000 | still has information
lag:  11 | p_value: 0.0001 | still has information
lag:  12 | p_value: 0.0001 | still has information
lag:  13 | p_value: 0.0001 | still has information
lag:  14 | p_value: 0.0001 | still has information
lag:  15 | p_value: 0.0000 | still has information
lag:  16 | p_value: 0.0001 | still has information
lag:  17 | p_value: 0.0001 | still has information
lag:  18 | p_value: 0.0001 | still has information
lag:  19 | p_value: 0.0002 | still has information
lag:  20 | p_value: 0.0003 | still has information
lag:  21 | p_value: 0.0002 | still has information
lag:  22 | p_value: 0.0003 | still has information
lag:  23 | p_value: 0.0005 | still has information
lag:  24 | p_value: 0.0004 | still has information
lag:  25 | p_value: 0.0005 | still has information
lag:  26 | p_value: 0.0008 | still has information
lag:  27 | p_value: 0.0005 | still has information
lag:  28 | p_value: 0.0007 | still has information
lag:  29 | p_value: 0.0004 | still has information
lag:  30 | p_value: 0.0004 | still has information
lag:  31 | p_value: 0.0002 | still has information
lag:  32 | p_value: 0.0003 | still has information
lag:  33 | p_value: 0.0000 | still has information
lag:  34 | p_value: 0.0000 | still has information
lag:  35 | p_value: 0.0000 | still has information
lag:  36 | p_value: 0.0000 | still has information
lag:  37 | p_value: 0.0000 | still has information
lag:  38 | p_value: 0.0000 | still has information
lag:  39 | p_value: 0.0000 | still has information
lag:  40 | p_value: 0.0000 | still has information

H = train['VIXCM30_Premium'].tolist()
ARIMA_ORDER = (3,0,4)
ARMA_prem = []
ARMA_PRED=[]
VIX = []
ACTUAL_FUTURE=[]
ACTUAL_prem = []
n = 0

for d,row in train.tail(1).iterrows():
    VIX +=[row['VIX']]

for d,row in test.iterrows():
    arima_model = ARIMA(H, order=ARIMA_ORDER)
    arima_fit = arima_model.fit()
    y_hat = arima_fit.predict(start=len(H), end=len(H) + 5)[0]
    ARMA_prem += [y_hat]
    ARMA_PRED+=[y_hat+VIX[-1]]
    VIX+=[row['VIX']]
    
    H += [row['VIXCM30_Premium']]
    ACTUAL_prem += [row['VIXCM30_Premium']]
    ACTUAL_FUTURE+=[row['VIXCM30']]
    n += 1


# plotting the first 200 points
plt.plot(ARMA_PRED[:200], label='predicted')
plt.plot(ACTUAL_FUTURE[:200], color='b', label='actual')
plt.legend()

<matplotlib.legend.Legend at 0x1c2e3019b0>

rolling_ar_errors2 = measure_error(ACTUAL_FUTURE, ARMA_PRED, label='Rolling ARIMA(%d,%d,%d)' % ARIMA_ORDER)

RMSE = 0.8361463719379066
RMSE_pcent = 0.04739375291785126
label = Rolling ARIMA(3,0,4)

ar_cm2 = plot_confusion_matrix(ACTUAL_FUTURE, ARMA_PRED)

Confusion Matrix
Raw
[[185 123]
 [165 115]]
Normalized
[[0.31462585 0.20918367]
 [0.28061224 0.19557823]]

make_money(ACTUAL_FUTURE, ARMA_PRED)

Total profit or loss: 15.6848
Number of correct days: 302
total incorrect days: 286

data	model	First lag in Ljung Box	RMSE_Percent	PnL
VIXCM30_Premium	ARMA(2,1)	5	0.04744	+\$ 1.56300
VIXCM30_Premium	ARMA(3,4)	10	0.04739	+\$ 15.6848

We can see that this order model also has information left in the lagged residuals (but starting at 10 rather than 5) and it also performs slightly better than the ARMA(2,1) model. This makes us a little cautious, but let's check if this problem persists in the ARMA models for the other Premium variables as well. We will spot check with the VIXCM60 Premium.

ARMA Model for VIXCM60 Premium¶

First let's take a look at the PACF and ACF:

X2= train['VIXCM60_Premium']

auto = acf(X2)
print(auto)
y = plot_acf(X2, lags=60)
plt.show()
y = plot_pacf(X2, lags=40)
plt.show()

[1.         0.95290316 0.92320006 0.90125888 0.88073219 0.86684725
 0.8488799  0.8317255  0.82183554 0.80894865 0.79688858 0.77657838
 0.75628295 0.74004291 0.72610523 0.70945886 0.69358428 0.67398217
 0.6540956  0.63959913 0.62170176 0.6032417  0.59016855 0.57442543
 0.56071146 0.54849787 0.53332192 0.52152519 0.50556092 0.48949031
 0.47956356 0.46535086 0.44882695 0.43283906 0.41861135 0.41088356
 0.40227172 0.39214396 0.38157287 0.36778869 0.35501799]

We will run a grid search to help us decide what p and q values to use in our model:

grid_result = run_grid_over_arma(X2, p_lags=list(range(5)), q_lags=range(5))

['p(0)', 'n/a', '11908.22', '10764.37', ' 9956.62', ' 9580.24']
['p(1)', ' 8299.45', ' 8222.75', ' 8208.13', ' 8213.54', ' 8213.48']
['p(2)', ' 8240.40', ' 8202.76', ' 8210.53', ' 8218.27', ' 8219.33']
['p(3)', ' 8220.75', ' 8210.53', ' 8217.39', ' 8224.26', ' 8211.52']
['p(4)', ' 8223.52', ' 8218.26', ' 8223.89', 'n/a', 'n/a']
+------+----------+----------+----------+----------+----------+
| p\q  |  q (0)   |  q (1)   |  q (2)   |  q (3)   |  q (4)   |
+------+----------+----------+----------+----------+----------+
| p(0) |   n/a    | 11908.22 | 10764.37 |  9956.62 |  9580.24 |
| p(1) |  8299.45 |  8222.75 |  8208.13 |  8213.54 |  8213.48 |
| p(2) |  8240.40 |  8202.76 |  8210.53 |  8218.27 |  8219.33 |
| p(3) |  8220.75 |  8210.53 |  8217.39 |  8224.26 |  8211.52 |
| p(4) |  8223.52 |  8218.26 |  8223.89 |   n/a    |   n/a    |
+------+----------+----------+----------+----------+----------+
bad values: [(0, 0), (4, 3), (4, 4)]
time required 11.767310857772827

b = sorted(list(grid_result.items()), key=lambda k: k[1])
print(b)

[((2, 1), 8202.764414557952), ((1, 2), 8208.130009167802), ((2, 2), 8210.527132456586), ((3, 1), 8210.527162769662), ((3, 4), 8211.523492331788), ((1, 4), 8213.482916654084), ((1, 3), 8213.537119955572), ((3, 2), 8217.394445200434), ((4, 1), 8218.258294896628), ((2, 3), 8218.268074030846), ((2, 4), 8219.32718617729), ((3, 0), 8220.754381057428), ((1, 1), 8222.75177438938), ((4, 0), 8223.519608762306), ((4, 2), 8223.885593759196), ((3, 3), 8224.257522342912), ((2, 0), 8240.398729208286), ((1, 0), 8299.445989273678), ((0, 4), 9580.237456651997), ((0, 3), 9956.617795085402), ((0, 2), 10764.365937165772), ((0, 1), 11908.215104282828)]

Again we will try an ARMA(2,1) model:

model = ARIMA(X2, order=(2,0,1))
r = model.fit()
r.summary()
print('RMSE (in sample):',np.sqrt(np.mean(r.sigma2)))

RMSE (in sample): 1.3711005884788794

train2=train['VIXCM60_Premium'].dropna().tolist()
arima_model2 = ARIMA(train2, order=(2,0,1))
model_fit2 = arima_model2.fit()

We can check to see if there is any information left in the residuals:

run_ljung_box(model_fit2.resid,lags=40)

Ljung-Box Test
lag:   7 | p_value: 0.0113 | still has information
lag:   8 | p_value: 0.0072 | still has information
lag:   9 | p_value: 0.0108 | still has information
lag:  10 | p_value: 0.0000 | still has information
lag:  11 | p_value: 0.0000 | still has information
lag:  12 | p_value: 0.0000 | still has information
lag:  13 | p_value: 0.0000 | still has information
lag:  14 | p_value: 0.0000 | still has information
lag:  15 | p_value: 0.0000 | still has information
lag:  16 | p_value: 0.0000 | still has information
lag:  17 | p_value: 0.0000 | still has information
lag:  18 | p_value: 0.0000 | still has information
lag:  19 | p_value: 0.0000 | still has information
lag:  20 | p_value: 0.0000 | still has information
lag:  21 | p_value: 0.0000 | still has information
lag:  22 | p_value: 0.0000 | still has information
lag:  23 | p_value: 0.0000 | still has information
lag:  24 | p_value: 0.0000 | still has information
lag:  25 | p_value: 0.0000 | still has information
lag:  26 | p_value: 0.0000 | still has information
lag:  27 | p_value: 0.0000 | still has information
lag:  28 | p_value: 0.0000 | still has information
lag:  29 | p_value: 0.0000 | still has information
lag:  30 | p_value: 0.0000 | still has information
lag:  31 | p_value: 0.0000 | still has information
lag:  32 | p_value: 0.0000 | still has information
lag:  33 | p_value: 0.0000 | still has information
lag:  34 | p_value: 0.0000 | still has information
lag:  35 | p_value: 0.0000 | still has information
lag:  36 | p_value: 0.0000 | still has information
lag:  37 | p_value: 0.0000 | still has information
lag:  38 | p_value: 0.0000 | still has information
lag:  39 | p_value: 0.0000 | still has information
lag:  40 | p_value: 0.0000 | still has information

Unfortunately, there is still a lot of data left in the residuals. But let's just check to see how the ARMA performs out of sample regardless.

H = train['VIXCM60_Premium'].tolist()
ARIMA_ORDER = (2,0,1)
ARMA_prem2 = []
ARMA_PRED2=[]
VIX = []
ACTUAL_FUTURE2=[]
ACTUAL_prem2 = []
n = 0

for d,row in train.tail(1).iterrows():
    VIX +=[row['VIX']]

for d,row in test.iterrows():
    arima_model = ARIMA(H, order=ARIMA_ORDER)
    arima_fit = arima_model.fit()
    y_hat = arima_fit.predict(start=len(H), end=len(H) + 5)[0]
    ARMA_prem2 += [y_hat]
    ARMA_PRED2+=[y_hat+VIX[-1]]
    VIX+=[row['VIX']]
    
    H += [row['VIXCM60_Premium']]
    ACTUAL_prem2 += [row['VIXCM60_Premium']]
    ACTUAL_FUTURE2+=[row['VIXCM60']]
    n += 1


# plotting the first 200 points
plt.plot(ARMA_PRED2[:200], label='predicted')
plt.plot(ACTUAL_FUTURE2[:200], color='b', label='actual')
plt.legend()

<matplotlib.legend.Legend at 0x1c2d5ad5c0>

rolling_ar_errors = measure_error(ACTUAL_FUTURE2, ARMA_PRED2, label='Rolling ARIMA(%d,%d,%d)' % ARIMA_ORDER)

RMSE = 0.6519801144622868
RMSE_pcent = 0.03534142555932097
label = Rolling ARIMA(2,0,1)

ar_cm = plot_confusion_matrix(ACTUAL_FUTURE2, ARMA_PRED2)

Confusion Matrix
Raw
[[185 116]
 [184 103]]
Normalized
[[0.31462585 0.19727891]
 [0.31292517 0.17517007]]

make_money(ACTUAL_FUTURE2, ARMA_PRED2)

Total profit or loss: -2.8920
Number of correct days: 288
total incorrect days: 300

Now let's try our ARMA (3,4) model to see if we can solve the issue of the information that is remaining in the residuals:

model2 = ARIMA(X, order=(3,0,4))
r2 = model2.fit()
r2.summary()
print('RMSE (in sample):',np.sqrt(np.mean(r2.sigma2)))

RMSE (in sample): 1.1897402525654754

arima_model = ARIMA(train1, order=(3,0,4))
model_fit2 = arima_model.fit()

run_ljung_box(r2.resid,lags=40)

Ljung-Box Test
lag:  10 | p_value: 0.0000 | still has information
lag:  11 | p_value: 0.0001 | still has information
lag:  12 | p_value: 0.0001 | still has information
lag:  13 | p_value: 0.0001 | still has information
lag:  14 | p_value: 0.0001 | still has information
lag:  15 | p_value: 0.0000 | still has information
lag:  16 | p_value: 0.0001 | still has information
lag:  17 | p_value: 0.0001 | still has information
lag:  18 | p_value: 0.0001 | still has information
lag:  19 | p_value: 0.0002 | still has information
lag:  20 | p_value: 0.0003 | still has information
lag:  21 | p_value: 0.0002 | still has information
lag:  22 | p_value: 0.0003 | still has information
lag:  23 | p_value: 0.0005 | still has information
lag:  24 | p_value: 0.0004 | still has information
lag:  25 | p_value: 0.0005 | still has information
lag:  26 | p_value: 0.0008 | still has information
lag:  27 | p_value: 0.0005 | still has information
lag:  28 | p_value: 0.0007 | still has information
lag:  29 | p_value: 0.0004 | still has information
lag:  30 | p_value: 0.0004 | still has information
lag:  31 | p_value: 0.0002 | still has information
lag:  32 | p_value: 0.0003 | still has information
lag:  33 | p_value: 0.0000 | still has information
lag:  34 | p_value: 0.0000 | still has information
lag:  35 | p_value: 0.0000 | still has information
lag:  36 | p_value: 0.0000 | still has information
lag:  37 | p_value: 0.0000 | still has information
lag:  38 | p_value: 0.0000 | still has information
lag:  39 | p_value: 0.0000 | still has information
lag:  40 | p_value: 0.0000 | still has information

H = train['VIXCM60_Premium'].tolist()
ARIMA_ORDER = (3,0,4)
ARMA_prem = []
ARMA_PRED=[]
VIX = []
ACTUAL_FUTURE=[]
ACTUAL_prem = []
n = 0

for d,row in train.tail(1).iterrows():
    VIX +=[row['VIX']]

for d,row in test.iterrows():
    arima_model = ARIMA(H, order=ARIMA_ORDER)
    arima_fit = arima_model.fit()
    y_hat = arima_fit.predict(start=len(H), end=len(H) + 5)[0]
    ARMA_prem += [y_hat]
    ARMA_PRED+=[y_hat+VIX[-1]]
    VIX+=[row['VIX']]
    
    H += [row['VIXCM60_Premium']]
    ACTUAL_prem += [row['VIXCM60_Premium']]
    ACTUAL_FUTURE+=[row['VIXCM60']]
    n += 1


# plotting the first 200 points
plt.plot(ARMA_PRED[:200], label='predicted')
plt.plot(ACTUAL_FUTURE[:200], color='b', label='actual')
plt.legend()

<matplotlib.legend.Legend at 0x1c3449e668>

rolling_ar_errors2 = measure_error(ACTUAL_FUTURE, ARMA_PRED, label='Rolling ARIMA(%d,%d,%d)' % ARIMA_ORDER)

RMSE = 0.6562020202035773
RMSE_pcent = 0.03557027942183079
label = Rolling ARIMA(3,0,4)

ar_cm2 = plot_confusion_matrix(ACTUAL_FUTURE, ARMA_PRED)

Confusion Matrix
Raw
[[190 111]
 [184 103]]
Normalized
[[0.32312925 0.18877551]
 [0.31292517 0.17517007]]

make_money(ACTUAL_FUTURE, ARMA_PRED)

Total profit or loss: -6.8196
Number of correct days: 293
total incorrect days: 295

data	model	First lag in Ljung Box	RMSE_Percent	PnL
VIXCM60_Premium	ARMA(2,1)	7	0.03534	-\$ 2.8920
VIXCM60_Premium	ARMA(3,4)	10	0.03557	-\$ 6.8196

Again we see that when using the ARMA model, we find there is still information left in the residuals. This makes us concerned that our ARMA models are not properly capturing all of the information available, but we feel a bit more confident that our AR models can. Because of this, and the fact that our AR models look to be making us more money, we will choose to use the AR models to develop our strategy, rather than continue running more ARMA models (which is great since we prefer simpler models and ARMA models take much longer to run).

Revised Strategy
¶

Revising Our Buy/Sell Strategy¶

We decided to check how our PnL might change if we allow the sale or purchase of any given future to take over the 5 next days, rather than requiring this action to happen the very next day. To do this, we created a new function, shown below, where the reccommendation (buy +1 or sell -1) follows the same strategy as before, buy if VIXCM30 is predicted to be higher tomorrow, sell if it is predicted to be lower. Then we will multiply the direction by $Future Price_{t*}$ - $Future Price_{t}$. Where $t*$ is the optimal day (max or min) over the next 5 trading days.

For example, if the VIXCM30 was 20 today and the predicted VIXCM30 over the next 5 days were 21,22,24,19,20 respectively, our strategy would tells us to buy (+1) since 20<21 and then we would sell in 3 days because this is the date with the predicted highest VIXCM30 price (24). Then let's day the real price on that 3rd say ends up being 23. Our PnL would be (+1)*(23-20) = +3.

def make_money_week(actual,predicted):
    d = {'Predicted': predicted,
     'Actual': actual}
    PnL_df = pd.DataFrame(d)
    
    Total_PnL =0
    direction=0
    correct=0
    incorrect=0
    
    for i in range(0, (len(PnL_df)-1)):
        if PnL_df['Predicted'][i+1][0]>PnL_df['Actual'][i][0]:
            direction=1
            best_day=max(PnL_df['Predicted'][i+1])
            pos=np.unravel_index(np.argmax((PnL_df['Predicted'][i+1]), axis=None), (PnL_df['Predicted'][i+1]).shape)
            actual_sell=PnL_df['Actual'][i][pos]
        else:
            direction=-1
            best_day=min(PnL_df['Predicted'][i+1])
            pos=np.unravel_index(np.argmin((PnL_df['Predicted'][i+1]), axis=None), (PnL_df['Predicted'][i+1]).shape)
            actual_sell=PnL_df['Actual'][i][pos]
        daily_PnL= direction*(actual_sell-PnL_df['Actual'][i][0])
        Total_PnL+=daily_PnL
        if daily_PnL>=0:
            correct+=1
        else:
            incorrect+=1
            
    print('Total profit or loss: %.4f' % (Total_PnL))
    print('Number of correct days: %d' % (correct))
    print('total incorrect days: %d' % (incorrect))

ACTUAL_LIST_week=[]
for i in range(0,(len(test)-4)):
    new=[test.VIXCM30.iloc[i],test.VIXCM30.iloc[i+1],
         test.VIXCM30.iloc[i+2],test.VIXCM30.iloc[i+3],test.VIXCM30.iloc[i+4]]
    new=np.asarray(new)
    ACTUAL_LIST_week+=[new]

H = train['VIXCM30_Premium'].tolist()
AR_ORDER = 49
AR_PRED_week_lists=[]
VIX=[]
for d,row in train.tail(1).iterrows():
    VIX+=[row['VIX']]
for d,row in test.iterrows():
    ar_model = AR(H)
    ar_fit = ar_model.fit(maxlag=AR_ORDER)
    y_hat = ar_fit.predict(start=len(H), end=len(H) + 4)
    pred_list=[]
    pred=y_hat+VIX[-1]
    pred_list.append(pred)
    AR_PRED_week_lists+=pred_list
    VIX+=[row['VIX']] 
    H += [row['VIXCM30_Premium']]

make_money_week(ACTUAL_LIST_week,AR_PRED_week_lists[:-4])

Total profit or loss: 20.5517
Number of correct days: 336
total incorrect days: 248

ACTUAL_LIST_week60=[]
for i in range(0,(len(test)-4)):
    new=[test.VIXCM60.iloc[i],test.VIXCM60.iloc[i+1],
         test.VIXCM60.iloc[i+2],test.VIXCM60.iloc[i+3],test.VIXCM60.iloc[i+4]]
    new=np.asarray(new)
    ACTUAL_LIST_week60+=[new]

H = train['VIXCM60_Premium'].tolist()
AR_ORDER = 35
AR_PRED_week_lists60=[]
VIX=[]
for d,row in train.tail(1).iterrows():
    VIX+=[row['VIX']]
for d,row in test.iterrows():
    ar_model = AR(H)
    ar_fit = ar_model.fit(maxlag=AR_ORDER)
    y_hat = ar_fit.predict(start=len(H), end=len(H) + 4)
    pred_list=[]
    pred=y_hat+VIX[-1]
    pred_list.append(pred)
    AR_PRED_week_lists60+=pred_list
    VIX+=[row['VIX']] 
    H += [row['VIXCM60_Premium']]

make_money_week(ACTUAL_LIST_week60,AR_PRED_week_lists60[:-4])

Total profit or loss: 8.9814
Number of correct days: 349
total incorrect days: 235

ACTUAL_LIST_week90=[]
for i in range(0,(len(test)-4)):
    new=[test.VIXCM90.iloc[i],test.VIXCM90.iloc[i+1],
         test.VIXCM90.iloc[i+2],test.VIXCM90.iloc[i+3],test.VIXCM90.iloc[i+4]]
    new=np.asarray(new)
    ACTUAL_LIST_week90+=[new]

H = train['VIXCM90_Premium'].tolist()
AR_ORDER = 35
AR_PRED_week_lists90=[]
VIX=[]
for d,row in train.tail(1).iterrows():
    VIX+=[row['VIX']]
for d,row in test.iterrows():
    ar_model = AR(H)
    ar_fit = ar_model.fit(maxlag=AR_ORDER)
    y_hat = ar_fit.predict(start=len(H), end=len(H) + 4)
    pred_list=[]
    pred=y_hat+VIX[-1]
    pred_list.append(pred)
    AR_PRED_week_lists90+=pred_list
    VIX+=[row['VIX']] 
    H += [row['VIXCM90_Premium']]

make_money_week(ACTUAL_LIST_week90,AR_PRED_week_lists90[:-4])

Total profit or loss: 12.4234
Number of correct days: 348
total incorrect days: 236

ACTUAL_LIST_week120=[]
for i in range(0,(len(test)-4)):
    new=[test.VIXCM120.iloc[i],test.VIXCM120.iloc[i+1],
         test.VIXCM120.iloc[i+2],test.VIXCM120.iloc[i+3],test.VIXCM120.iloc[i+4]]
    new=np.asarray(new)
    ACTUAL_LIST_week120+=[new]

H = train['VIXCM120_Premium'].tolist()
AR_ORDER = 35
AR_PRED_week_lists120=[]
VIX=[]
for d,row in train.tail(1).iterrows():
    VIX+=[row['VIX']]
for d,row in test.iterrows():
    ar_model = AR(H)
    ar_fit = ar_model.fit(maxlag=AR_ORDER)
    y_hat = ar_fit.predict(start=len(H), end=len(H) + 4)
    pred_list=[]
    pred=y_hat+VIX[-1]
    pred_list.append(pred)
    AR_PRED_week_lists120+=pred_list
    VIX+=[row['VIX']] 
    H += [row['VIXCM120_Premium']]

make_money_week(ACTUAL_LIST_week120,AR_PRED_week_lists120[:-4])

Total profit or loss: 4.5694
Number of correct days: 351
total incorrect days: 233

ACTUAL_LIST_week150=[]
for i in range(0,(len(test)-4)):
    new=[test.VIXCM150.iloc[i],test.VIXCM150.iloc[i+1],
         test.VIXCM150.iloc[i+2],test.VIXCM150.iloc[i+3],test.VIXCM150.iloc[i+4]]
    new=np.asarray(new)
    ACTUAL_LIST_week150+=[new]

H = train['VIXCM150_Premium'].tolist()
AR_ORDER = 49
AR_PRED_week_lists150=[]
VIX=[]
for d,row in train.tail(1).iterrows():
    VIX+=[row['VIX']]
for d,row in test.iterrows():
    ar_model = AR(H)
    ar_fit = ar_model.fit(maxlag=AR_ORDER)
    y_hat = ar_fit.predict(start=len(H), end=len(H) + 4)
    pred_list=[]
    pred=y_hat+VIX[-1]
    pred_list.append(pred)
    AR_PRED_week_lists150+=pred_list
    VIX+=[row['VIX']] 
    H += [row['VIXCM150_Premium']]

make_money_week(ACTUAL_LIST_week150,AR_PRED_week_lists150[:-4])

Total profit or loss: 10.3753
Number of correct days: 363
total incorrect days: 221

data	model	PnL (over a week)
VIXCM30_Premium	AR(49)	+\$ 20.5517
VIXCM60_Premium	AR(35)	+\$ 8.9814
VIXCM90_Premium	AR(35)	+\$ 12.4234
VIXCM120_Premium	AR(35)	+\$ 4.5694
VIXCM150_Premium	AR(49)	+\$ 10.3753

Using this new strategy, we end up with higher profits for each model!

Conclusion
¶

After trying ARCH, GARCH, VECM, AR, and ARMA models we came to find that our AR models resulted in the most profitable strategy. It is worth noting that the RSMEs for ARCH and GARCH models were the lowest, so they were the best fit, but these models were predicting the market volitality rather than the price of the futures. In regards to predicting the VIX futures' prices, AR models did the best, but they still didn't result in a strong edge. When revising our strategy to allow the sale or purchase of a future over the 5 next days, rather than requiring this action to happen the very next day, we improved our profits.

	Date	SPY	RealVol	VIX	VIXCM30	VIXCM60	VIXCM90	VIXCM120	VIXCM150
0	3/23/2005	92.578594	9.120305	14.06	14.506154	15.065909	15.454545	15.739375	15.955000
1	3/24/2005	92.689358	6.457921	13.42	14.225263	14.894773	15.242500	15.524375	15.768125
2	3/28/2005	92.823878	5.578653	13.75	14.398649	14.993636	15.320909	15.586250	15.820625
3	3/29/2005	92.206709	9.290422	14.49	14.931667	15.219773	15.458409	15.667813	15.864688
4	3/30/2005	93.512316	7.313855	13.64	14.428571	15.091364	15.398182	15.600312	15.764375

	SPY	RealVol	VIX	VIXCM30	VIXCM60	VIXCM90	VIXCM120	VIXCM150
count	2943.000000	2943.000000	2943.000000	2943.000000	2943.000000	2943.000000	2915.000000	509.000000
mean	132.905896	11.228967	19.586290	20.737151	21.419516	21.851063	22.277360	17.815324
std	41.630454	8.782778	9.433398	7.751411	6.942745	6.461391	6.149277	4.051598
min	58.030194	1.725931	9.890000	11.763684	12.900000	13.569474	14.047656	14.298125
25%	102.141627	6.420053	13.455000	15.270175	16.291053	16.963577	17.662237	15.432500
50%	118.316563	8.551970	16.670000	18.522917	19.500000	20.193421	20.660526	16.255079
75%	170.156606	12.769802	22.535000	24.030250	24.680550	25.076158	25.684722	18.905000
max	221.520000	110.821685	80.860000	62.706364	54.465500	48.518947	45.061053	37.655000

	SPY	RealVol	VIX	VIXCM30	VIXCM60	VIXCM90	VIXCM120	VIXCM150
Date
3/23/2005	92.578594	9.120305	14.06	14.506154	15.065909	15.454545	15.739375	15.955000
3/24/2005	92.689358	6.457921	13.42	14.225263	14.894773	15.242500	15.524375	15.768125
3/28/2005	92.823878	5.578653	13.75	14.398649	14.993636	15.320909	15.586250	15.820625
3/29/2005	92.206709	9.290422	14.49	14.931667	15.219773	15.458409	15.667813	15.864688
3/30/2005	93.512316	7.313855	13.64	14.428571	15.091364	15.398182	15.600312	15.764375

	coef	std err	z	P>\|z\|	[0.025	0.975]
exog1	-6.127e-05	0.000	-0.151	0.880	-0.001	0.001
L1.VIX	-0.3103	0.042	-7.325	0.000	-0.393	-0.227
L1.VIXCM30	0.2802	0.088	3.169	0.002	0.107	0.453
L2.VIX	-0.2080	0.043	-4.863	0.000	-0.292	-0.124
L2.VIXCM30	0.1658	0.089	1.872	0.061	-0.008	0.339
L3.VIX	-0.1166	0.043	-2.718	0.007	-0.201	-0.033
L3.VIXCM30	0.0563	0.089	0.635	0.525	-0.117	0.230
L4.VIX	-0.1284	0.043	-2.989	0.003	-0.213	-0.044
L4.VIXCM30	0.0175	0.089	0.197	0.843	-0.156	0.191
L5.VIX	0.0449	0.043	1.045	0.296	-0.039	0.129
L5.VIXCM30	-0.2715	0.089	-3.066	0.002	-0.445	-0.098
L6.VIX	-0.0783	0.043	-1.825	0.068	-0.162	0.006
L6.VIXCM30	0.1179	0.089	1.330	0.183	-0.056	0.292
L7.VIX	-0.0661	0.043	-1.543	0.123	-0.150	0.018
L7.VIXCM30	-0.0227	0.089	-0.256	0.798	-0.196	0.151
L8.VIX	0.0640	0.043	1.495	0.135	-0.020	0.148
L8.VIXCM30	-0.1825	0.089	-2.062	0.039	-0.356	-0.009
L9.VIX	0.0832	0.043	1.944	0.052	-0.001	0.167
L9.VIXCM30	-0.2264	0.089	-2.555	0.011	-0.400	-0.053
L10.VIX	0.1476	0.043	3.450	0.001	0.064	0.231
L10.VIXCM30	-0.1111	0.089	-1.252	0.211	-0.285	0.063
L11.VIX	0.1314	0.043	3.065	0.002	0.047	0.215
L11.VIXCM30	-0.2368	0.089	-2.670	0.008	-0.411	-0.063
L12.VIX	-0.0393	0.043	-0.918	0.359	-0.123	0.045
L12.VIXCM30	0.1105	0.089	1.245	0.213	-0.063	0.284
L13.VIX	-0.0114	0.043	-0.267	0.789	-0.095	0.072
L13.VIXCM30	-0.0027	0.089	-0.030	0.976	-0.177	0.171
L14.VIX	0.1181	0.042	2.787	0.005	0.035	0.201
L14.VIXCM30	-0.2756	0.089	-3.112	0.002	-0.449	-0.102
L15.VIX	0.0855	0.042	2.023	0.043	0.003	0.168
L15.VIXCM30	-0.1531	0.089	-1.730	0.084	-0.327	0.020
L16.VIX	0.1011	0.042	2.393	0.017	0.018	0.184
L16.VIXCM30	-0.0786	0.088	-0.888	0.375	-0.252	0.095
L17.VIX	0.0745	0.042	1.772	0.076	-0.008	0.157
L17.VIXCM30	-0.0325	0.088	-0.368	0.713	-0.206	0.141
L18.VIX	0.0057	0.042	0.136	0.892	-0.076	0.087
L18.VIXCM30	-0.0607	0.088	-0.689	0.491	-0.233	0.112
L19.VIX	0.1579	0.041	3.845	0.000	0.077	0.238
L19.VIXCM30	-0.3513	0.088	-4.007	0.000	-0.523	-0.179
L20.VIX	0.0813	0.040	2.047	0.041	0.003	0.159
L20.VIXCM30	-0.1109	0.087	-1.278	0.201	-0.281	0.059

	coef	std err	z	P>\|z\|	[0.025	0.975]
beta.1	1.0000	0	0	0.000	1.000	1.000
beta.2	-1.0485	0.024	-44.311	0.000	-1.095	-1.002
lin_trend	0.0003	0.000	0.605	0.545	-0.001	0.001

	SPY	RealVol	VIX	VIXCM30	VIXCM60	VIXCM90	VIXCM120	VIXCM150	spy_ret	variance
Date
2005-03-23	92.578594	9.120305	14.06	14.506154	15.065909	15.454545	15.739375	15.955000	NaN	NaN
2005-03-24	92.689358	6.457921	13.42	14.225263	14.894773	15.242500	15.524375	15.768125	0.001196	6.802115e-07
2005-03-28	92.823878	5.578653	13.75	14.398649	14.993636	15.320909	15.586250	15.820625	0.001451	1.165579e-06
2005-03-29	92.206709	9.290422	14.49	14.931667	15.219773	15.458409	15.667813	15.864688	-0.006649	4.928736e-05
2005-03-30	93.512316	7.313855	13.64	14.428571	15.091364	15.398182	15.600312	15.764375	0.014160	1.901057e-04

	SPY	RealVol	VIX	VIXCM30	VIXCM60	VIXCM90	VIXCM120	VIXCM150	spy_ret	variance	VIXCM30_Premium	VIXCM60_Premium	VIXCM90_Premium	VIXCM120_Premium	VIXCM150_Premium
Date
2014-07-23	189.716544	3.895231	11.52	13.471053	14.396	15.210526	15.990909	16.771292	0.002220	3.416250e-06	1.951053	2.876	3.690526	4.470909	5.251292
2014-07-24	189.726090	3.964469	11.84	13.692105	14.578	15.310526	16.072727	16.834928	0.000050	1.032740e-07	1.852105	2.738	3.470526	4.232727	4.994928
2014-07-25	188.837875	5.437378	12.69	14.005263	14.806	15.476316	16.211364	16.946411	-0.004682	2.553529e-05	1.315263	2.116	2.786316	3.521364	4.256411
2014-07-28	188.914283	6.551183	12.56	13.897368	14.726	15.368421	16.136364	16.904306	0.000405	1.085231e-09	1.337368	2.166	2.808421	3.576364	4.344306
2014-07-29	188.102462	6.450938	13.28	14.002632	14.700	15.331579	16.118182	16.904785	-0.004297	2.179938e-05	0.722632	1.420	2.051579	2.838182	3.624785

Det. terms outside the coint. relation & lagged endog. parameters for equation VIXCM30
	coef	std err	z	P>\|z\|	[0.025	0.975]
exog1	0.0005	0.000	2.424	0.015	8.68e-05	0.001
L1.VIX	-0.0867	0.020	-4.446	0.000	-0.125	-0.048
L1.VIXCM30	0.0775	0.041	1.904	0.057	-0.002	0.157
L2.VIX	-0.0532	0.020	-2.704	0.007	-0.092	-0.015
L2.VIXCM30	0.0270	0.041	0.662	0.508	-0.053	0.107
L3.VIX	-0.0202	0.020	-1.022	0.307	-0.059	0.019
L3.VIXCM30	-0.0199	0.041	-0.488	0.626	-0.100	0.060
L4.VIX	-0.0497	0.020	-2.511	0.012	-0.088	-0.011
L4.VIXCM30	0.0118	0.041	0.290	0.772	-0.068	0.092
L5.VIX	0.0002	0.020	0.012	0.990	-0.039	0.039
L5.VIXCM30	-0.1161	0.041	-2.849	0.004	-0.196	-0.036
L6.VIX	-0.0098	0.020	-0.498	0.618	-0.049	0.029
L6.VIXCM30	0.0058	0.041	0.143	0.886	-0.074	0.086
L7.VIX	-0.0090	0.020	-0.454	0.650	-0.048	0.030
L7.VIXCM30	-0.0396	0.041	-0.970	0.332	-0.120	0.040
L8.VIX	0.0070	0.020	0.357	0.721	-0.032	0.046
L8.VIXCM30	-0.0480	0.041	-1.179	0.239	-0.128	0.032
L9.VIX	0.0188	0.020	0.952	0.341	-0.020	0.057
L9.VIXCM30	-0.0691	0.041	-1.694	0.090	-0.149	0.011
L10.VIX	0.0602	0.020	3.058	0.002	0.022	0.099
L10.VIXCM30	-0.0520	0.041	-1.273	0.203	-0.132	0.028
L11.VIX	0.0457	0.020	2.314	0.021	0.007	0.084
L11.VIXCM30	-0.0869	0.041	-2.127	0.033	-0.167	-0.007
L12.VIX	-0.0147	0.020	-0.744	0.457	-0.053	0.024
L12.VIXCM30	0.0169	0.041	0.415	0.678	-0.063	0.097
L13.VIX	-0.0245	0.020	-1.246	0.213	-0.063	0.014
L13.VIXCM30	0.0353	0.041	0.863	0.388	-0.045	0.115
L14.VIX	0.0320	0.020	1.639	0.101	-0.006	0.070
L14.VIXCM30	-0.1189	0.041	-2.916	0.004	-0.199	-0.039
L15.VIX	0.0244	0.019	1.255	0.209	-0.014	0.063
L15.VIXCM30	-0.0405	0.041	-0.993	0.321	-0.120	0.039
L16.VIX	0.0622	0.019	3.199	0.001	0.024	0.100
L16.VIXCM30	-0.0757	0.041	-1.859	0.063	-0.156	0.004
L17.VIX	0.0355	0.019	1.835	0.067	-0.002	0.073
L17.VIXCM30	-0.0265	0.041	-0.652	0.515	-0.106	0.053
L18.VIX	0.0274	0.019	1.427	0.154	-0.010	0.065
L18.VIXCM30	-0.0672	0.041	-1.657	0.097	-0.147	0.012
L19.VIX	0.0711	0.019	3.759	0.000	0.034	0.108
L19.VIXCM30	-0.1493	0.040	-3.698	0.000	-0.228	-0.070

Det. terms outside the coint. relation & lagged endog. parameters for equation VIX
	coef	std err	z	P>\|z\|	[0.025	0.975]
exog1	-0.0002	0.000	-0.358	0.720	-0.001	0.001
L1.VIX	-0.3074	0.038	-8.014	0.000	-0.383	-0.232
L1.VIXCM60	0.3885	0.111	3.510	0.000	0.172	0.605
L2.VIX	-0.1803	0.039	-4.649	0.000	-0.256	-0.104
L2.VIXCM60	0.1506	0.111	1.357	0.175	-0.067	0.368
L3.VIX	-0.0317	0.039	-0.816	0.414	-0.108	0.044
L3.VIXCM60	-0.1912	0.111	-1.722	0.085	-0.409	0.026
L4.VIX	-0.0966	0.039	-2.486	0.013	-0.173	-0.020
L4.VIXCM60	-0.0289	0.111	-0.261	0.794	-0.246	0.188
L5.VIX	0.0369	0.039	0.950	0.342	-0.039	0.113
L5.VIXCM60	-0.3169	0.111	-2.861	0.004	-0.534	-0.100
L6.VIX	-0.0321	0.039	-0.825	0.409	-0.108	0.044
L6.VIXCM60	0.0223	0.111	0.201	0.841	-0.195	0.240
L7.VIX	-0.0779	0.039	-2.004	0.045	-0.154	-0.002
L7.VIXCM60	0.0392	0.111	0.354	0.724	-0.178	0.257
L8.VIX	0.0485	0.039	1.251	0.211	-0.028	0.125
L8.VIXCM60	-0.1901	0.111	-1.715	0.086	-0.407	0.027
L9.VIX	0.0451	0.039	1.163	0.245	-0.031	0.121
L9.VIXCM60	-0.1629	0.111	-1.469	0.142	-0.380	0.054
L10.VIX	0.1432	0.039	3.687	0.000	0.067	0.219
L10.VIXCM60	-0.1457	0.111	-1.314	0.189	-0.363	0.072
L11.VIX	0.0310	0.039	0.795	0.427	-0.045	0.107
L11.VIXCM60	-0.0177	0.111	-0.160	0.873	-0.235	0.200
L12.VIX	-0.0805	0.039	-2.070	0.038	-0.157	-0.004
L12.VIXCM60	0.2636	0.111	2.380	0.017	0.047	0.481
L13.VIX	-0.0061	0.039	-0.156	0.876	-0.082	0.070
L13.VIXCM60	-0.0395	0.111	-0.356	0.722	-0.257	0.178
L14.VIX	0.1011	0.039	2.616	0.009	0.025	0.177
L14.VIXCM60	-0.3335	0.111	-3.013	0.003	-0.550	-0.117
L15.VIX	0.0588	0.039	1.524	0.127	-0.017	0.135
L15.VIXCM60	-0.1454	0.111	-1.312	0.189	-0.363	0.072
L16.VIX	0.1134	0.039	2.942	0.003	0.038	0.189
L16.VIXCM60	-0.1718	0.111	-1.553	0.121	-0.389	0.045
L17.VIX	0.0410	0.038	1.067	0.286	-0.034	0.116
L17.VIXCM60	0.0449	0.111	0.406	0.685	-0.172	0.262
L18.VIX	-0.0064	0.038	-0.168	0.867	-0.081	0.069
L18.VIXCM60	-0.0661	0.111	-0.598	0.550	-0.283	0.151
L19.VIX	0.1312	0.038	3.471	0.001	0.057	0.205
L19.VIXCM60	-0.3803	0.110	-3.446	0.001	-0.597	-0.164

Det. terms outside the coint. relation & lagged endog. parameters for equation VIXCM60
	coef	std err	z	P>\|z\|	[0.025	0.975]
exog1	0.0003	0.000	2.344	0.019	5.55e-05	0.001
L1.VIX	-0.0794	0.013	-6.094	0.000	-0.105	-0.054
L1.VIXCM60	0.1398	0.038	3.722	0.000	0.066	0.213
L2.VIX	-0.0521	0.013	-3.958	0.000	-0.078	-0.026
L2.VIXCM60	0.0745	0.038	1.977	0.048	0.001	0.148
L3.VIX	-0.0023	0.013	-0.174	0.862	-0.028	0.024
L3.VIXCM60	-0.0637	0.038	-1.691	0.091	-0.138	0.010
L4.VIX	-0.0309	0.013	-2.343	0.019	-0.057	-0.005
L4.VIXCM60	0.0140	0.038	0.372	0.710	-0.060	0.088
L5.VIX	-0.0048	0.013	-0.361	0.718	-0.031	0.021
L5.VIXCM60	-0.0818	0.038	-2.176	0.030	-0.156	-0.008
L6.VIX	0.0054	0.013	0.407	0.684	-0.020	0.031
L6.VIXCM60	-0.0359	0.038	-0.954	0.340	-0.110	0.038
L7.VIX	-0.0083	0.013	-0.633	0.527	-0.034	0.018
L7.VIXCM60	-0.0190	0.038	-0.505	0.613	-0.093	0.055
L8.VIX	0.0031	0.013	0.236	0.813	-0.023	0.029
L8.VIXCM60	-0.0286	0.038	-0.760	0.447	-0.102	0.045
L9.VIX	0.0059	0.013	0.449	0.653	-0.020	0.032
L9.VIXCM60	-0.0260	0.038	-0.690	0.490	-0.100	0.048
L10.VIX	0.0527	0.013	3.999	0.000	0.027	0.079
L10.VIXCM60	-0.0686	0.038	-1.821	0.069	-0.142	0.005
L11.VIX	0.0095	0.013	0.719	0.472	-0.016	0.035
L11.VIXCM60	-0.0117	0.038	-0.310	0.757	-0.085	0.062
L12.VIX	-0.0376	0.013	-2.845	0.004	-0.063	-0.012
L12.VIXCM60	0.0905	0.038	2.406	0.016	0.017	0.164
L13.VIX	-0.0093	0.013	-0.706	0.480	-0.035	0.017
L13.VIXCM60	0.0183	0.038	0.487	0.626	-0.055	0.092
L14.VIX	0.0246	0.013	1.873	0.061	-0.001	0.050
L14.VIXCM60	-0.1169	0.038	-3.110	0.002	-0.191	-0.043
L15.VIX	0.0174	0.013	1.324	0.185	-0.008	0.043
L15.VIXCM60	-0.0452	0.038	-1.201	0.230	-0.119	0.029
L16.VIX	0.0538	0.013	4.111	0.000	0.028	0.079
L16.VIXCM60	-0.0924	0.038	-2.460	0.014	-0.166	-0.019
L17.VIX	0.0284	0.013	2.176	0.030	0.003	0.054
L17.VIXCM60	-0.0115	0.038	-0.306	0.759	-0.085	0.062
L18.VIX	0.0143	0.013	1.100	0.271	-0.011	0.040
L18.VIXCM60	-0.0613	0.038	-1.633	0.102	-0.135	0.012
L19.VIX	0.0482	0.013	3.757	0.000	0.023	0.073
L19.VIXCM60	-0.1233	0.037	-3.292	0.001	-0.197	-0.050

Det. terms outside the coint. relation & lagged endog. parameters for equation VIX
	coef	std err	z	P>\|z\|	[0.025	0.975]
exog1	-0.0003	0.000	-0.687	0.492	-0.001	0.001
L1.VIX	-0.3300	0.036	-9.224	0.000	-0.400	-0.260
L1.VIXCM90	0.6012	0.123	4.897	0.000	0.361	0.842
L2.VIX	-0.1274	0.036	-3.499	0.000	-0.199	-0.056
L2.VIXCM90	-0.0456	0.123	-0.370	0.711	-0.287	0.196
L3.VIX	-0.0130	0.036	-0.356	0.722	-0.084	0.059
L3.VIXCM90	-0.2710	0.123	-2.197	0.028	-0.513	-0.029
L4.VIX	-0.0855	0.036	-2.346	0.019	-0.157	-0.014
L4.VIXCM90	-0.0476	0.123	-0.386	0.700	-0.289	0.194
L5.VIX	0.0202	0.036	0.553	0.580	-0.051	0.092
L5.VIXCM90	-0.2812	0.123	-2.282	0.023	-0.523	-0.040
L6.VIX	-0.0125	0.037	-0.341	0.733	-0.084	0.059
L6.VIXCM90	-0.0312	0.123	-0.253	0.800	-0.273	0.211
L7.VIX	-0.0775	0.037	-2.120	0.034	-0.149	-0.006
L7.VIXCM90	0.0171	0.124	0.138	0.890	-0.225	0.259
L8.VIX	0.0228	0.037	0.623	0.533	-0.049	0.094
L8.VIXCM90	-0.1166	0.123	-0.945	0.345	-0.359	0.125
L9.VIX	0.0387	0.037	1.059	0.290	-0.033	0.110
L9.VIXCM90	-0.1416	0.123	-1.147	0.251	-0.383	0.100
L10.VIX	0.1620	0.037	4.434	0.000	0.090	0.234
L10.VIXCM90	-0.2312	0.124	-1.871	0.061	-0.473	0.011
L11.VIX	0.0099	0.037	0.269	0.788	-0.062	0.082
L11.VIXCM90	0.0697	0.123	0.565	0.572	-0.172	0.312
L12.VIX	-0.0693	0.037	-1.890	0.059	-0.141	0.003
L12.VIXCM90	0.2686	0.123	2.180	0.029	0.027	0.510
L13.VIX	0.0069	0.037	0.188	0.851	-0.065	0.079
L13.VIXCM90	-0.0892	0.123	-0.724	0.469	-0.331	0.153
L14.VIX	0.1127	0.036	3.099	0.002	0.041	0.184
L14.VIXCM90	-0.4459	0.123	-3.623	0.000	-0.687	-0.205
L15.VIX	0.0487	0.036	1.340	0.180	-0.023	0.120
L15.VIXCM90	-0.1032	0.123	-0.839	0.401	-0.344	0.138
L16.VIX	0.1249	0.036	3.450	0.001	0.054	0.196
L16.VIXCM90	-0.2495	0.123	-2.035	0.042	-0.490	-0.009
L17.VIX	0.0165	0.036	0.457	0.648	-0.054	0.087
L17.VIXCM90	0.1683	0.123	1.372	0.170	-0.072	0.409
L18.VIX	-0.0113	0.036	-0.314	0.754	-0.082	0.059
L18.VIXCM90	-0.0790	0.123	-0.645	0.519	-0.319	0.161
L19.VIX	0.0603	0.036	1.696	0.090	-0.009	0.130
L19.VIXCM90	-0.1572	0.123	-1.283	0.200	-0.397	0.083

VIX Project : Forecasting Market Volatility¶

Table of Contents¶

Data notes:¶

Project goals:¶

First we will import the data and define functions which will be used throughout the notebook:¶

Importing the data¶

Before we jump in we will make sure we understand the data available:¶

Now we will explore the data:¶

Addressing NAs¶

Before We Build Models We Will Establish How We Will Measure Their Success:¶

Models To Be Tested:¶

ARCH and GARCH Models¶

Time to check if we made money with this strategy:¶

Summary of ARCH and GARCH Models¶

VECM Models¶

VIXCM30¶

VIXCM60¶

VIXCM90¶

VIXCM120¶

VIXCM150¶

Summary of VECM Models:¶

AR Models¶

Defining New Premium Variables:¶

Now we will check if all of the premium variables are stationary:¶

AR Model for VIXCM30 Premium¶

Check to see how much money could be made out of sample:¶

AR Model for VIXCM60 Premium¶

Check to see how much money could be made out of sample:¶

AR Model for VIXCM90 Premium¶

Check to see how much money could be made out of sample:¶

VIXCM120 Premium¶

Check to see how much money could be made out of sample:¶

VIXCM150 Premium¶

Check to see how much money could be made out of sample:¶

Summary of AR Models:¶

ARMA Models¶

ARMA Model for VIXCM30 Premium¶

ARMA Model for VIXCM60 Premium¶

Revising Our Buy/Sell Strategy¶

Table of Contents
¶

VECM Models
¶

AR Models
¶